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Abstract 

Due  to  the  versatility  of  its  strueture,  the  semi-Markov  proeess  is  a  powerful  modeling 
tool  used  to  deseribe  eomplex  systems.  Though  similar  in  structure  to  continuous  time 
Markov  chains,  semi-Markov  processes  allow  for  any  transition  time  distribution  which 
enables  these  processes  to  fit  a  wider  range  of  problems  than  the  continuous  time  Markov 
chain.  While  semi-Markov  processes  have  been  applied  in  fields  as  varied  as  biostatistics 
and  finance,  there  does  not  exist  a  theoretically-based,  systematic  method  to  determine  if  a 
semi-Markov  process  accurately  fits  the  underlying  data  used  to  create  the  model. 

In  fields  such  as  regression  and  analysis  of  variance,  the  quality  of  the  predictive 
model  is  judged  in  part  by  the  goodness  of  fit  of  the  model  which  relates  the  expected 
observation  values  with  the  actual  observations.  A  similar  methodology  for  semi-Markov 
processes  would  provide  immediate  insight  in  the  efficacy  of  the  fitted  model  and  would 
allow  competing  models  to  be  directly  compared  with  one  another. 

This  dissertation  presents  a  methodology  to  measure  the  adequacy  of  a  fitted  semi- 
Markov  process.  To  this  end,  a  technique  to  assess  the  likelihood  that  a  data  sample  could 
be  generated  by  a  specific  semi-Markov  process  is  developed,  including  a  newly  proposed 
goodness  of  fit  metric.  This  technique  relies  on  the  covariance  structure  of  the  semi-Markov 
process;  thus,  a  method  to  estimate  the  covariance  structure  is  also  proposed.  The  technique 
is  applied  to  real  and  simulated  data  to  demonstrate  the  goodness  of  fit  metric’s  utility  in 
model  validation  and  its  ability  to  identify  potential  covariate  factors  within  the  model. 
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TESTING  THE  ADEQUACY  OE  A  SEMI-MARKOV  PROCESS 


I.  Introduction 

The  area  of  stoehastie  proeesses  known  as  multi-state  modeling  deseribes  systems  in 
a  vast  array  of  fields,  from  finanee  and  health  care  to  game  theory  and  reliability.  Typically, 
these  models  decompose  a  complex  system  into  a  series  of  connected  states  based  on 
observed  data.  Although  observed  data  provides  the  basis  for  the  model,  a  modeler’s 
judgment  also  plays  a  substantial  role  in  model  development.  The  variability  inherent  in 
observational  data  and  the  modeler’s  biases  may  result  in  a  potentially  inadequate  system 
model.  Other  modeling  arenas,  such  as  regression  and  distribution  fitting,  utilize  various 
quality  metrics  to  confirm  model  adequacy;  an  analogous  technique  does  not  exist  for 
general  multi-state  modeling. 

Simple  quality  metrics  determine  the  likelihood  that  an  observed  sample  comes  from 
a  specific  model,  as  illustrated  by  the  analysis  of  variance  E-statistic  and  Pearson’s  Chi- 
squared  goodness  of  fit  statistic.  More  advanced  metrics,  such  as  R^,  assess  the  degree 
of  fit  between  the  model  and  the  observed  data.  All  of  these  metrics  are  based  in  theory, 
simple  to  compute,  and  provide  the  modeler  with  immediate  insight  regarding  the  validity 
of  a  model. 

Currently,  multi- state  modeling  efforts  omit  any  reference  to  model  goodness  of  fit  [1- 
3].  This  is  because  the  computational  simplicity  is  believed  to  be  lost  due  to  the  complexity 
of  the  model.  In  reality  though,  certain  forms  of  quality  metrics  may  be  applied  to  multi¬ 
state  models.  However,  past  efforts  have  seen  limited  success  [4]  due  to  the  subjective 
nature  of  many  of  the  available  metrics  and  nuances  in  data  collection.  Additionally, 
quality  assessment  tools  often  fail  to  provide  definitive  results  in  multi-state  models  due 
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to  data  collection  techniques.  Specifically,  censoring  within  data  collection  methods  means 
that  actual  state  transition  times  are  not  observed  and  must  be  approximated,  which  adds 
additional  uncertainty  in  any  quality  metric  calculations.  These  drawbacks  result  in  quality 
metrics  losing  their  allure,  especially  when  the  modeler  must  test  multiple  competing 
models.  To  be  of  use  for  multi-state  models,  the  quality  metric  must  be  able  to  differentiate 
between  models  in  a  quick  and  reliable  manner. 

The  work  of  this  dissertation  utilizes  a  Chi-squared  test  statistic  based  on  a  model’s 
covariance  structure  to  create  a  goodness  of  fit  metric  for  a  semi-Markov  process.  The 
metric  is  developed  and  applied  to  a  semi-Markov  process  since  the  semi-Markov  process 
is  a  generalized  class  that  contains  any  Markov  process  including  discrete-time  and 
continuous-time  Markov  chains.  This  metric  is  used  to  test  the  adequacy  of  the  proposed 
model  with  respect  to  the  observed  data  sample.  The  technique  to  calculate  the  metric 
and  test  model  adequacy  incorporates  a  variety  of  nuances  contained  within  the  structure 
of  a  semi-Markov  process  in  order  to  elicit  the  required  expectations  for  the  metric.  The 
testing  technique  is  demonstrated  through  a  series  of  well-defined  models.  Additionally, 
an  application  of  this  technique  demonstrates  its  usefulness  in  model  selection  by  deciding 
whether  or  not  potential  covariate  factors  significantly  impact  the  model’s  performance. 

This  dissertation  is  structured  into  seven  chapters.  Following  this  first  introductory 
chapter,  background  information  is  provided  in  Chapter  2  to  describe  the  multi-state  model 
including  the  semi-Markov  process.  Chapter  3  provides  the  theory  for  calculating  the 
covariance  of  the  Markov  renewal  process  which  is  a  surrogate  for  the  covariance  of  a  semi- 
Markov  process.  Chapter  4  presents  the  goodness  of  fit  metric  and  its  properties  as  well  as 
provides  a  simulation  study  to  demonstrate  the  robustness  of  the  proposed  metric.  Chapter 
5  demonstrates  how  the  goodness  of  fit  metric  may  be  used  to  test  the  model  structure  for 
the  presence  of  a  possible  covariate.  Finally,  Chapter  6  demonstrates  the  application  of 
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the  goodness  of  fit  metric  to  examine  trends  in  the  tempo  of  Body  Mass  Index  values  in 
childhood  and  conclusions  are  provided  in  Chapter  7. 
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II.  Background 


2.1  Stochastic  Processes 

A  stochastic  process  is  a  collection  of  random  variables  indexed  on  some  set  T .  The 
evolution  of  these  random  variables  ean  be  used  to  represent  a  system  ehanging  over  time. 
For  example,  let  X{t)  be  the  priee  of  a  barrel  of  oil  at  time  t,  then  X  =  {X(t),  t  e  T}  is  a 
stoehastie  process  that  deseribes  the  price  evolution  of  a  barrel  of  oil  over  time  T .  T  is 
also  known  as  the  index  set.  This  index  set  ean  be  diserete  (the  elosing  priee  eaeh  day)  or 
eontinuous  (the  instanteous  priee  throughout  the  day).  An  aetual  realization  of  X  is  known 
as  a  sample  path. 

Stoehastie  processes  trace  the  evolution  of  a  system  as  it  assumes  various  states,  5.  A 
state,  say  j'l  or  S2,  is  an  observable  loeation  of  the  stoehastie  proeess  and  the  state  space 
(S  =  {j'l,  S2, ...})  is  the  full  set  of  possible  states.  Although  S  eould  be  M,  in  this  work 
S  will  be  limited  to  a  diserete  state  spaee.  Figure  2.1  shows  a  four  state  system  with 
potential  transitions  between  states  indieated  by  arrows.  Notice  the  system  never  returns 
after  departing  State  1  and  never  leaves  after  arriving  at  State  4. 


Figure  2.1:  Basic  state  model 
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2.2  Markov  Processes 


A  Markov  process  is  a  stochastic  process  for  which  the  Markov  property  holds  [5]. 
The  Markov  property  is,  given  the  entire  history  of  the  process,  the  conditional  distribution 
of  a  future  state  is  only  dependent  on  the  current  state  [6].  The  Markov  property  reduces 
the  required  system  memory  to  just  the  current  state  rather  than  an  entire  history  of  the 
system.  Instead  of  tracking  every  state  visited  and  the  time  spent  in  every  state,  a  Markov 
process  only  needs  to  track  the  current  system  state.  From  a  computational  standpoint,  this 
represents  a  significant  reduction  in  physical  memory  requirements. 

Markov  processes  may  take  on  many  forms  depending  on  the  structure  of  the  system. 
A  Markov  process  with  a  discrete  state  space  is  known  as  a  Markov  chain.  A  Markov  chain 
is  fully  determined  by  its  matrix  of  transition  probabilities  from  State  i  directly  to  State  j 
for  each  i  and  y  in  S .  With  this  transition  probability  matrix,  P,  and  the  current  state,  many 
features  of  the  system  can  be  determined,  including  the  likelihood  that  the  system  will  be  in 
a  particular  state  after  the  next  transition  or  after  the  n-th  transition,  how  many  transitions 
will  occur  until  a  specific  state  is  reached,  and  the  amount  of  time  until  the  system  returns 
to  the  current  state.  Although  the  state  space  is  discrete,  the  index  set  for  a  Markov  chain 
can  be  either  discrete  or  continuous  based  upon  how  time  is  measured  and  tracked. 

This  leads  to  the  definition  of  a  discrete  time  Markov  chain.  Let  denote  the  system 
state  after  the  n-th  transition. 

Definition  Discrete  Time  Markov  Chain  [7]  A  stochastic  process  {X„,n  >  0}  with 
countable  state  space  5  is  a  discrete  time  Markov  chain  if  the  following  hold  for  all  integer¬ 
valued  n>  0,  and  for  each  i,  j  in  5 : 

1.  Xn  G  S,  and 

2.  =  i\Xn  =  kXn-U  ...,Xo)  =  =  i\Xn  =  i). 

Item  2  is  the  Markov  property  with  respect  to  the  discrete  time  Markov  chain.  Note  that 
time  is  not  explicitly  tracked  in  the  discrete  time  Markov  chain.  In  this  type  of  process. 
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each  transition  is  considered  a  single  time  step.  In  effect,  the  n  is  the  number  of  time  steps 
and  can  be  related  back  to  the  actual  index  set,  T,  for  the  process. 

Definition  Continuous  Time  Markov  Chain  [7]  A  stochastic  process  >  0}  with 

countable  state  space  S  is  a  continuous  time  Markov  chain  if  the  following  are  satisfied  by 
the  sequence  {Aq,  (A„,  ¥„),  n>  1}  for  all  integer- valued  n  >  0,  and  for  each  i,  j  inS: 

P(Xn+l  =  j,Yn+\  >  y\Xn  =  i,Yfi,Xfi-l,Yfi-\,  ■..,Xi,Yi,Xq) 

=  P{X,=j,Y,  >y\X^  =  i)  =  pije-^^y. 

In  this  formulation,  A„  is  the  system  state  after  the  n-th  transition,  is  the  time  the  system 
requires  to  transition  from  A„_i  to  A„,  also  known  as  the  sojourn  time,  pij  is  the  probability 
of  transitioning  directly  from  State  i  to  State  j,  and  qi  is  the  exponential  parameter  for  the 
sojourn  time  distribution  for  State  i  such  that  qik  where  k  is  any  state  in  S . 

In  a  continuous  time  Markov  chain,  every  transition  sojourn  time  follows  an 
Exponential  distribution  which  maintains  the  Markov  property  due  to  the  memoryless 
nature  of  the  distribution.  Further,  the  transition  probability  matrix,  P,  is  equal  to  [p,;]. 
The  typical  continuous  time  Markov  chain  is  formulated  so  that  pu  =  0,  for  all  i.  This 
formulation  eliminates  immediate  transitions  back  to  the  current  state  which  simplifies 
moment  calculations.  If  a  system  does  allow  self  transitions,  additional  states  can  be  added 
in  order  to  represent  the  self  transitions  without  altering  the  performance  of  the  system. 
Often,  the  ^,y’s  are  placed  in  a  matrix  Q  =  [^,y]  which  is  known  as  the  generator  matrix, 
where  qu  =  -  Qik  =  -qi-  As  a  result,  the  row  sums  in  Q  are  equal  to  0.  A  continuous 
time  Markov  chain  is  fully  described  by  Q. 

The  continuous  time  Markov  chain  construct  is  restricted  to  modeling  sojourn  times 
using  exponential  distributions.  To  allow  any  sojourn  time  distribution,  more  complex 
models,  such  as  a  semi-Markov  process,  must  be  used.  The  semi-Markov  process  was 
originally  defined  by  Levy  [8],  Takacs  [9],  and  Smith  [10]  and  further  refined  by  Pyke  [11]. 
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Definition  Semi-Markov  Process  [7]  A  stochastic  process  >  0}  with  countable 

state  space  5  is  a  semi-Markov  process  if  the  sequence  {Aq,  (X„,  Yn),n  >  1}  satisfies  the 
following  for  all  integer-valued  n  >  0,  and  for  each  i,  j  in  S: 

P(Xn+l  =  j,Yn+\  <y\Xn  =  i,Yn,Xn-l,Yn-\,--.,Xi,Yi,XQ) 

=  P{X,=  pY,<y\X^  =  i)  =  Gij{y). 

Here,  A„  and  have  the  same  definition  as  a  continuous  time  Markov  chain,  and 

Gij  =  PijFij<y)^ 

which  is  the  probability  of  transitioning  from  State  i  directly  to  State  j  (pij)  times  the 
cumulative  distribution  function  of  the  sojourn  time  distribution  between  States  i  and  j 
(Fijiy)).  The  semi-Markov  process  kernel  matrix  is  given  by  G(y)  =  [G;y(y)]  for  i,  j  in 
S,y>0. 

As  the  formal  definition  implies,  the  semi-Markov  process  is  similar  to  the  continuous 
time  Markov  chain.  The  semi-Markov  process  consists  of  a  transition  probability  matrix 
and  distributions  for  the  sojourn  time  between  transitions.  However,  by  relaxing  the 
exponential  sojourn  distribution  requirement,  the  Markov  property  no  longer  applies  at 
every  time  t,  but  only  at  the  actual  transition  times  [5].  For  example,  suppose  a  system 
transitions  into  State  1  at  time  t  and  will  transition  into  State  2  at  some  time  t  -v  s.  If 
the  system  was  a  continuous  time  Markov  chain,  the  transition  time  would  follow  an 
Exponential(/l)  distribution  and  for  any  time  t  -l-  /,  i  <  s,  the  expected  transition  time  would 
be  A.  On  the  other  hand,  if  the  system  was  a  semi-Markov  process,  the  transition  time 
would  follow  some  general  distribution  that  may  not  adhere  to  the  memoryless  property, 
and  the  expected  transition  time  would  not  be  the  same  for  any  time  t  -v  i,  i  <  s.  Therefore, 
the  semi-Markov  process  must  track  the  history  between  transitions  by  maintaining  the 
time  since  the  last  transition. 

The  semi-Markov  process  construct  is  useful  in  providing  added  fidelity  to  real  world 
systems  because  of  the  flexibility  in  choosing  sojourn  distributions.  Although  semi-Markov 
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processes  are  commonly  associated  with  reliability  and  survival  analyses  [12],  they  have 
also  been  applied  to  a  variety  of  problems  in  areas  such  as  insurance  [2],  biology  [13], 
finance  [1],  and  computer  science  [14].  The  basic  approach  to  modeling  with  a  semi- 
Markov  process  from  observed  data  is  the  following: 

1 .  Define  the  state  space,  S . 

2.  Estimate  the  transition  probability  distribution  matrix,  P. 

3.  Estimate  the  cumulative  distribution  functions  of  the  sojourn  times  for  F. 

As  a  semi-Markov  process  evolves  over  time,  the  number  of  visits  to  each  state  forms 
a  Markov  renewal  process  as  shown  in  Pyke  [11]. 

Definition  Markov  Renewal  Process  [7]  A  process  N  is  called  a  Markov  renewal  process, 
if  the  vector  N  =  {N{t)  =  [A/0]»  ^  ^  0,  for  all  j  in  S}, 

where  Nj(t)  records  the  number  of  visits  the  system  has  made  to  a  State  j  at  time  t  given  an 
initial  state  at  time  0.  Eigure  2.2  illustrates  the  relationship  between  semi-Markov  processes 
and  Markov  renewal  processes  by  showing  the  transitions  for  the  semi-Markov  process 
on  the  left  and  the  counting  process  for  State  3  on  the  right.  States  1,  2,  and  4  have 
similar  Markov  renewal  process  graphs.  While  the  semi-Markov  process  tracks  a  single 
state  occupancy  value,  i.e.  the  system  is  currently  in  State  4,  the  Markov  renewal  process 
tracks  the  number  of  visits  to  each  of  the  S  states.  The  link  between  the  semi-Markov 
process  and  the  Markov  renewal  process  is  exploited  in  the  creation  of  the  goodness  of  fit 
metric  in  Chapter  4  as  the  metric  calculations  are  performed  on  the  Markov  renewal  process 
instead  of  the  actual  semi-Markov  process. 

The  Markov  renewal  process  has  a  well-defined  expectation.  As  illustrated  in  Pyke 
[15],  the  expected  count  in  the  Eaplace  domain  at  a  fixed  for  each  state  is 

X(E[A(0])  =  (I-q(0r'q(0,  (2.1) 
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Figure  2.2:  Semi-Markov  proeess  and  related  Markov  renewal  proeess  for  State  3 


where  q(0  =  P  o  f(t),  o  is  the  Hadamard  product  of  two  matrices,  and  f  is  a  matrix  of 
probability  density  functions  that  correspond  to  F.  For  an  observed  semi-Markov  process, 
the  corresponding  Markov  renewal  process,  N,  has  an  expectation,  A  =E[Aj],  for  all  jinS . 
The  vector  value  M  =  N  -  A  represents  the  unexplained  variation  or  error  of  the  Markov 
renewal  process  which  forms  a  Martingale. 

A  martingale,  M(t),  is  a  stochastic  process  such  that  E[A/(t)]  <  oo  for  all  t  and 
E[M(t)\M{s)]  =  M{s)  for  s  <  t.  Random  walks  and  classical  Brownian  motion  are  two 
common  examples  of  martingales.  Submartingales  and  supermartingales  are  two  extended 
classes  of  martingales.  They  are  both  stochastic  processes  with  E[M(t)]  <  oo,  V  t, 
but  a  submartingale  has  E[M{t)\M{s)]  >  M{s),  s  <  t,  while  a  supermartingale  has 
E[M(t)\M(s)]  <  M{s),  s  <  t.  A  true  martingale  is,  in  fact,  a  special  case  of  both  a 
submartingale  and  a  supermartingale.  Martingales,  submartingales,  and  supermartingales 
all  belong  to  a  larger  class  called  semimartingales  [16]. 

2.3  Semi-Markov  Process  Diagnostics 

In  most  modeling  realms,  diagnostic  tools  are  used  to  assess  the  model  fit  with 
respect  to  observed  data.  Eor  multi-state  models,  including  semi-Markov  processes,  model 
diagnostics  are  underdeveloped  when  compared  with  other  areas  such  as  linear  modeling. 
While  a  variety  of  diagnostic  tools  have  been  developed,  the  majority  only  work  with  very 
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specific  types  of  data  and  do  not  provide  a  clear  metric  for  whether  or  not  the  model  truly  fits 
the  data.  Titman  and  Sharpies  [4]  provide  a  comprehensive  review  of  the  various  techniques 
for  diagnosing  model  fit  in  multi-state  models.  Their  review  includes  the  shortcomings  of 
most  diagnostic  tools  with  respect  to  panel-observed  medical  data. 

Of  the  diagnostic  tools  reviewed,  prevalence  counts  and  Chi-squared  goodness  of  fit 
metrics  show  promise  in  assessing  model  fit  but  do  not  satisfy  the  requirements  of  a  formal 
diagnostic  test.  Prevalence  counts  [17]  compare  the  number  of  times  a  state  is  observed 
with  the  number  of  expected  observations  of  the  state  for  a  series  of  fixed  time.  The 
comparison  uses  the  metric 


Ht  = 


Et 


(2.2) 


where  O  is  the  observed  count  and  E  is  the  expected  count.  Large  values  of  H  indicate  poor 
model  fit,  but  prevalence  counts  cannot  formally  assess  model  fit  because  the  distribution 
of  H  is  unknown.  Under  a  balanced  observation  scheme,  meaning  the  observations  occur 
at  regular,  set  intervals. 


H 


ZEE 

t  i=l  7=1 


(Otij  -  Etij) 
Etij 


(2.3) 


follows  a  Chi-squared  distribution  where  t  is  the  set  of  observation  times,  i  is  the  occupied 
state  at  t,  and  j  is  the  occupied  state  at  t  -i-  1  [18].  Unfortunately,  balanced  observation 
schemes  are  rarely  attained  in  data  collection.  While  irregularly  observed  data  can  be  tested 
with  Equation  2.3  [19],  the  resulting  statistic  does  not  consistently  follow  the  Chi-squared 
distribution  in  all  cases. 

The  two  metrics  in  Equations  2.2  and  2.3  are  applications  of  the  generalized  Pearson 
Chi-squared  test  often  used  for  assessing  goodness  of  fit  [20].  The  test  statistic  is  created 
by  grouping  the  data  into  a  series  of  n  bins,  then  computing  the  following  value: 


E 

(=1 


(Oi  -  Eif 


(2.4) 


where  O  is  the  number  of  observed  instances  in  bin  i  and  E  is  the  expected  number  of 
instances  in  bin  i.  The  statistic  follows  a  Chi-squared  distribution  with  n  -  p  degrees 
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of  freedom  where  p  is  the  number  of  estimated  parameters  plus  1.  Structural  equation 
modeling  uses  the  upper  diagonal  of  a  covariance  matrix  such  that 


^  — Y — 

<•=1  j=i 


(2.5) 


where  Otj  and  Eij  are  the  observed  and  expected  covariances  between  States  i  and  j 
respectively  [21].  This  statistic  follows  a  Chi-squared  distribution  with  ^n{n  -l-  1)  degrees 
of  freedom. 

The  lack  of  a  true  formal  test  for  multi-state  models,  including  semi-Markov 
processes,  motivates  the  development  of  the  goodness  of  fit  metric  in  following 
chapters.  This  metric  builds  off  the  previous  attempts  outlined  in  Equations  2.2-2. 5  while 
incorporating  the  inherent  structure  of  the  semi-Markov  process  and  its  related  Markov 
renewal  process.  The  next  chapter  develops  the  covariance  structure  for  a  Markov  renewal 
process  as  a  surrogate  for  the  semi-Markov  process  covariance. 
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III.  The  Covariance  of  a  Markov  Renewal  Process 


The  Markov  renewal  proeess  N(t),  defined  in  Seetion  2.2,  is  represented  by  a  veetor 
eontaining  a  eount  of  the  number  of  transitions  into  eaeh  of  the  n  states  at  time  t.  A 
transition  probability  matrix,  P,  and  a  sojourn  time  distribution  matrix,  F,  uniquely  define 
a  speeifie  Markov  renewal  proeess.  While  the  moments  of  the  Markov  renewal  process 
are  well  defined,  see  Appendix  A,  there  is  no  explicit  formulation  for  the  covariance 
between  the  state  counts.  This  covariance,  27, y,  describes  how  the  counts  of  States  i  and 
j  within  the  Markov  renewal  process  change  together,  or  co-vary,  and  can  be  expressed  as 
a  function  of  the  correlation  between  the  states  and  the  individual  variances  of  the  states. 
Traditionally,  the  covariance  is  used  in  goodness  of  fit  testing  and  will  be  required  for  the 
metric  developed  in  Chapter  4.  This  chapter  describes  how  to  compute  the  covariance  for 
a  Markov  renewal  process  and  illustrates  the  computation  using  a  notional  semi-Markov 
process. 

3.1  Calculating  the  Markov  Renewal  Process  Covariance 

Let  N(t)  be  a  column  vector  of  state  counts,  A(t)  be  a  column  vector  of  expected  state 
counts,  and  M{t)  be  a  column  vector  of  state  Martingales  at  time  t  as  defined  in  Chapter  2. 
The  covariance  between  state  counts  is  then  given  by: 

2:[iV(t)]  =  Cov[A(t),  A(0]  =  E[(A(t)  -  E[A(0])(7V(t)  -  E[A(t)])'] 

=  E[(A(o  -  m)m)  -  my]  =  Emtwm  o.d 


12 


M,{tf  Mi(0M2(0  ... 

M,{t)M2{t)  M2(tf  ...  M2(t)MM 

MiiDMnit)  M2{t)Mn(i)  ...  Mnitf 

E[Mi(0"]  E[Mi(0M2(0]  ...  E[Mi(0M„(0] 

E[Mi(0M2(0]  E[M2(0']  ...  E[M2(0M„(0] 

E[Mi(0M„(0]  E[M2(0M„(0]  ...  E[M„(02] 

Var[M(0]  E[Mi(0M2(0]  ...  E[Mi(0M„(0] 

E[Mi(0M2(0]  Var[A^2(0]  ...  E[M2(0M„(0] 

E[Mi(OM„a)]  E[M2(0M„(0]  ...  Var[A^„(0] 

Thus  the  expectations  of  the  product  of  two  random  variables,  E[M,(t)My(t)],  comprise 

the  off-diagonal  portions  of  the  Markov  renewal  process  covariance  matrix.  While  the 
variance  of  N{t)  has  a  solution  in  the  Eaplace  domain  (Appendix  A),  E[M;(0My(0]  does 
not  at  this  time. 

Although  Z(t)  cannot  currently  be  calculated  directly,  simulation  will  provide  a 
reasonable  approximation.  This  simulation  involves  collecting  N  iterations  of  a  semi- 
Markov  process  out  to  a  target  time  t  and  calculating  the  Markov  renewal  processes,  Niit) 
where  i  goes  from  1  to  N.  The  simulated  covariance  for  the  Markov  renewal  process  is  then 
computed  as: 

m  =  "  E[A(t)])(AK0  -  E[N(t)]f  (3.2) 

based  on  the  covariance  calculation  for  a  sample  [22].  Eor  a  given  semi-Markov  process 
that  does  not  allow  self  transitions,  E[iV(t)]  of  the  related  Markov  renewal  process  [11]  is  a 
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constant  based  on  t.  In  the  Laplace  domain, 

X(E[iV(0|Xo  =  /])  =  [(I  -  q)-^q]  . ,  (3.3) 

where  Xo  is  the  initial  system  state,  ^  =  P  o  f(t)  at  each  t,  and  f  is  a  matrix  of  probability 
density  functions  that  correspond  to  F  [15].  An  Euler  inversion  technique  returns  the  time 
domain  expected  value. 


E[N(t)\Xo  =  i]  = 


Re  (X(E[a  +  iu]))  cos(ut)  du  , 


where  a  is  a  point  such  that  X(E[a])  is  continuous  to  the  right  and  Re(X(E[a  +  iu]))  is  the 
real  portion  of  the  complex  X(E[a  +  iu])  [23].  This  integral  can  be  calculated  numerically 
using  a  Eourier- series  summation.  Thus  only  the  number  of  iterations,  N,  and  the  time,  t, 
must  be  selected  by  the  user  to  estimate  the  simulated  covariance.  The  lower  bounds  on 
both  of  these  parameters  depend  on  the  semi-Markov  process  being  modeled. 

Eirst,  consider  a  four-state  semi-Markov  process  with  a  probability  transition  matrix 

Pn  Pi2  Pn  Pi4 
P21  P22  P23  P24 

P  =  (3.5) 

P3l  P32  P33  P34 

P4\  P42  P43  P44 

and  a  sojourn  distribution  matrix  of 
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Eigure  3.1  illustrates  the  generic  semi-Markov  process  defined  by  P  and  F  from 
Equations  3.5  and  3.6.  To  compute  the  covariance  matrix  of  this  fully  connected  semi- 
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Markov  process,  simulate  from  the  model  as  follows: 


0 

0 

0  . 

.  0 

0 

0 

0  . 

.  0 

N  = 

and  S[A(0]  = 

0 

5x1 

0 

0  . 

.  0 

T  =  0  and  i  =  0; 

Set  t  to  the  desired  length  of  the  simulation; 

Set  N  to  the  desired  number  of  simulation  iterations; 

Calculate  E[A/^(t)]  using  Equation  3.4; 

while  i  <N  do 

Select  an  initial  state,  Sc,  based  on  the  starting  distribution  of  the  system  (i.e. 
the  probability  that  the  systems  begins  in  a  particular  start); 

while  r  <  t  do 

Draw  a  random  variable,  Xi  ~Uniform(0,l); 

Use  P  to  find  the  next  state,  based  on  Xi  and  Sc', 

Draw  a  random  variable,  X2  ~Uniform(0,l); 

Eind  the  transition  time,  r„,  to  5  „  by  making  a  random  draw  from 
based  on  X2, 

T  =  T  +  Tn, 
if  r  <=  t  then 
Sc  —  S  n. 

Ns  =  Ns  +  1; 
end 
end 

mit)]  =  2:[iV(0]  +  J^,(N  -  E[Nm(N  -  E[iV(0])'; 

end 

Algorithm  1:  Covariance  Simulation 
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At  this  point,  N  equals  the  total  number  of  times  eaeh  state  was  visited  during  the 


simulation  prior  to  t. 


3.2  The  Covariance  of  the  Sample  Problem 

The  Covarianee  Simulation  is  demonstrated  for  a  sample  problem  using  assumed 
values.  This  proeess  does  not  allow  instantaneous  transitions  baek  to  the  eurrent  state, 
State  1  eannot  transition  direetly  to  State  3,  and  State  2  eannot  transition  direetly  to  State 
4.  Let 


0  0.75  0  0.25 

0.5  0  0.5  0 


(3.7) 


0.25  0.45  0  0.3 
0.15  0.35  0.5  0 


and 


0 

r(2,2) 

0 

r(2,2) 

r(2, 1) 

0 

r(0.5,4) 

0 

r(0.5, 1) 

r(0.5,4) 

0 

r(2,2) 

r(0.5,4) 

r(2,2) 

r(2, 1) 

0 

16 


where 


Jr>A: 

——u^-^e-Mu.  (3.9) 

0  r(a)yS“ 

Running  N  =10,000  iterations  of  the  Covariance  Simulation  at  t  =1,000  time  units  results 
in  the  following  estimated  covariance: 

59.297  38.523  -9.579  -6.323 

38.523  75.407  26.323  -15.108 

i:[iV(1000)]  =  .  (3.10) 

-9.579  26.323  70.904  18.971 

-6.323  -15.108  18.971  37.058 

The  true  variance  of  the  Markov  renewal  process  can  be  calculated  using  the  equation  for 
the  Var*[A^(t)]  found  in  Appendix  A  and 

Var[A(t)]  = -^  I  Re(Var*[a  + /m])cos(m0  dM.  (3.11) 

^  Jo 


Based  on  the  semi-Markov  process  variances,  the  true  2’[A(1000)]  diagonal  values  are 
59.2,  75.5,  70.8,  37.1  respectively.  The  estimates  produced  by  the  Covariance  Simulation 
are  almost  identical  to  the  true  values. 

The  choice  of  N  will  affect  the  estimate  of  2’(t).  Increasing  the  number  of  iterations 
used  to  estimate  the  covariance  matrix  will  result  in  a  closer  approximation,  but  at  the 
cost  of  increased  computational  time.  To  illustrate  the  effect  of  the  number  of  iterations, 
the  covariance  estimate  of  the  sample  problem  is  computed  at  t=  1,000  time  units  using 
N  =100,  500,  1,000,  5,000,  10,000,  50,000,  and  100,000  iterations,  repeated  across 
R=  1,000  samples.  The  average  estimate  of  the  covariance  matrix  and  its  corresponding 
95%  confidence  interval  around  each  Zjj  value  was  computed.  The  results  for  the  various 
iteration  sizes  are  as  follows: 
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000)]  and  (95%  confidence  interval)  based  on  A^=100  iterations  = 


59.8(47.2, 74.1) 
38.7(26.8,50.7) 
-9.9(-20.7,1.3) 
-6.4(- 14.5, 1.2) 


38.7(26.8,50.7) 
76.0(60.1,94.3) 
26.6(14.3,39.7) 
-15.1(-24.0,  -6.5) 


-9.9(-20.7, 1.3) 
26.6(14.3,39.7) 
71.7(55.6,88.4) 
19.2(10.8,28.1) 


-6.4(-14.5,1.2) 

-15.1(-24.0,-6.5) 

19.2(10.8,28.1) 

37.3(29.1,46.0) 


(3.12) 


2'[A^(1, 000)]  and  (95%  confidence  interval)  based  on  N-500  iterations  = 


59.1(52.7,65.5) 

38.6(33.1,44.5) 

-9.6(- 14.4, -4.4) 

-6.4(-10.1,-3.2) 

38.6(33.1,44.5) 

75.6(67.8,83.9) 

26.2(20.5,31.8) 

-15.3(-19.7,-11.1) 

-9.6(- 14.4, -4.4) 

26.2(20.5,31.8) 

70.9(64.0,78.3) 

19.0(15.0,23.1) 

-6.4(-10.1,-3.2) 

-15.3(-19.7,-11.1) 

19.0(15.0,23.1) 

37.2(33.3,41.3) 

2'[A^(1, 000)]  and  (95%  confidence  interval) 

based  on  1,000  iterations  = 

59.2(54.9,63.8) 

38.6(34.6,42.5) 

-9.7(-13.1,-6.4) 

-6.4(-8.9,-4.0) 

38.6(34.6,42.5) 

75.6(70.2,81.1) 

26.2(22.4, 29.9) 

-15.3(-18.2,-12.4) 

-9.7(-13.1,-6.4) 

26.2(22.4, 29.9) 

70.9(65.8,76.0) 

19.0(16.1,21.7) 

-6.4(-8.9,  -4.0) 

-15.3(-18.2,-12.4) 

19.0(16.1,21.7) 

37.2(34.5,40.1) 

2'[A^(1, 000)]  and  (95%  confidence  interval) 

based  on  A^=5,000  iterations  = 

59.1(57.3,61.0) 

38.5(36.8,40.2) 

-9.7(-ll.l,-8.2) 

-6.4(-7.5,-5.3) 

38.5(36.8,40.2) 

75.5(73.2,78.0) 

26.2(24.4, 27.9) 

-15.3(-16.5,-14.0) 

-9.7(-ll.l,-8.2) 

26.2(24.4, 27.9) 

70.9(68.5,73.2) 

19.0(17.8,20.2) 

-6.4(-7.5,  -5.3) 

-15.3(-16.5,-14.0) 

19.0(17.8,20.2) 

37.2(35.9,38.4) 

2'[A^(1, 000)]  and  (95%  confidence  interval) 

based  on  10,000 

iterations  = 

59.1(57.9,60.4) 

38.5(37.3,39.7) 

-9.6(-10.7,-8.5) 

-6.4(-7.2,  -5.7) 

38.5(37.3,39.7) 

75.5(73.7,77.3) 

26.2(25.0,27.5) 

-15.3(-16.2,-14.4) 

-9.6(- 10.7, -8.5) 

26.2(25.0,27.5) 

70.9(69.3,72.5) 

18.9(18.1,19.8) 

-6.4(-7.2,-5.7) 

-15.3(-16.2,-14.4) 

18.9(18.1,19.8) 

37.2(36.3,38.0) 
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2'[A^(1, 000)]  and  (95%  confidence  interval)  based  on  A^=50,000  iterations  = 


59.2(58.6,59.7) 

38.5(38.0,39.1) 

-9.7(-10.1,-9.2) 

-6.4(-6.8,-6.1) 


38.5(38.0,39.1) 

75.5(74.7,76.3) 

26.2(25.6,26.8) 

-15.3(-15.7,-14.9) 


-9.7(-10.1,-9.2) 

26.2(25.6,26.8) 

70.9(70.1,71.6) 

18.9(18.5,19.3) 


-6.4(-6.8,  -6.1) 
-15.3(-15.7,-14.9) 
18.9(18.5,19.3) 
37.2(36.8,37.5) 


2'[A^(1, 000)]  and  (95%  confidence  interval)  based  on  100,000  iterations  = 


59.2(58.7,59.6) 
38.5(38.1,38.9) 
-9.6(- 10.0, -9.3) 
-6.4(-6.6,  -6.2) 


38.5(38.1,38.9) 

75.5(74.9,76.1) 

26.2(25.8,26.6) 

-15.3(-15.5,-15.0) 


-9.6(-10.0,-9.3) 

26.2(25.8,26.6) 

70.9(70.3,71.4) 

18.9(18.7,19.2) 


-6.4(-6.6,  -6.2) 
-15.3(-15.5,-15.0) 
18.9(18.7,19.2) 
37.2(36.9,37.4) 


From  the  covariance  estimates,  the  average  diagonal  values  mirror  the  true  variance 
values  with  500  iterations,  albeit  with  wide  95%  confidence  intervals  of  10  units  which 
indicates  that  a  single  set  of  500  interations  may  not  provide  a  good  estimate  of  2’[A^(f)]. 
Figure  3.2  shows  the  decreasing  confidence  interval  sizes  as  the  number  of  iterations 
increase.  With  50,000  iterations,  the  confidence  intervals  drop  to  1  unit.  In  this  scenario, 
it  takes  the  same  amount  of  time  to  estimate  the  covariance  with  a  single  set  of  50,000 
iterations  as  it  does  to  average  1,000  samples  of  500  iterations  each.  However,  the  later 
method  will  provide  a  more  accurate  estimate  of  the  true  covariance  matrix  since  it  relies 
on  averaging  a  series  of  estimates  rather  than  using  in  single  point  estimate. 

While  t  is  dictated  by  the  observed  data  sample,  the  minimum  semi-Markov  process 
runtime  required  for  accurate  covariance  estimates  depends  largely  on  the  particular  model 
and  how  quickly  it  reaches  a  steady  state  on  average.  For  instance,  a  model  that  averages  10 
time  units  before  the  first  transition  cannot  have  an  accurate  covariance  estimate  for  N{5) 
because  the  Markov  renewal  process  counts  will  be  highly  dependent  on  the  initial  states. 
Meanwhile  a  small  model  that  averages  0.1  time  units  between  transitions  can  have  an 
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Figure  3.2:  State-to-State  95%  eonfidence  interval  widths  for  model  eovarianee 


an  aeeurate  eovarianee  estimate  for  N(5)  beeause  the  system  will  have  enough  transitions 
to  overeome  the  effeet  of  the  initial  state.  To  illustrate,  the  eorrelation  matrix  is  used  to 
normalize  the  eovarianee  estimates  across  time.  The  correlation  matrix  is  calculated  by 

Corr,y[iV(0]  =  2  (3-13) 

yizummzjjmt)] 
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Matrices  of  the  average  correlation  with  95%  confidence  intervals  around  each  Corr,y 
value  for  N(l),  N(20),  N(50),  and  N(IOO)  follow: 

Corr[A^(l)]  and  (95%  confidence  interval)  = 

1(0.89,1.1)  -0.05(-0.1 1,0.01)  0.02(-0.03,0.08)  0(-0.05, 0.06) 

-0.05(-0.1 1,0.01)  1(0.85,1.15)  0.33(0.24,0.44)  -0.05(-0.09, -0.02) 

0.02(-0.03,0.08)  0.33(0.24,0.44)  1(0.86,1.14)  0.03(-0.01,0.09) 

0(-0.05,0.06)  -0.05(-0.09,-0.02)  0.03(-0.01,0.09)  1(0.6,1.44) 

Corr[A^(20)]  and  (95%  confidence  interval)  = 

1(0.93,1.08)  0.5(0.44,0.56)  -0.12(-0.17, -0.07)  -0.12(-0.17, -0.07) 

0.5(0.44,0.56)  1(0.92,1.08)  0.37(0.31,0.44)  -0.28(-0.33, -0.23) 

-0.12(-0.17,-0.07)  0.37(0.31,0.44)  1(0.93,1.08)  0.32(0.27,0.38) 

-0.12(-0.17,-0.07)  -0.28(-0.33,-0.23)  0.32(0.27,0.38)  1(0.93,1.08) 

Corr[iV(50)]  and  (95%  confidence  interval)  = 

1(0.93,1.07)  0.55(0.49,0.61)  -0.14(-0.19, -0.09)  -0.13(-0.18, -0.08) 

0.55(0.49,0.61)  1(0.93,1.08)  0.36(0.31,0.42)  -0.28(-0.34, -0.23) 

-0.14(-0.19,-0.09)  0.36(0.31,0.42)  1(0.93,1.08)  0.35(0.3,0.41) 

-0.13(-0.18,-0.08)  -0.28(-0.34,-0.23)  0.35(0.3,0.41)  1(0.93,1.08) 

Corr[M100)]  and  (95%  confidence  interval)  = 

1(0.93,1.07)  0.56(0.5,0.63)  -0.14(-0.2, -0.09)  -0.13(-0.19,-0.08) 

0.56(0.5,0.63)  1(0.93,1.08)  0.36(0.3,0.42)  -0.29(-0.34, -0.24) 

-0.14(-0.2,-0.09)  0.36(0.3,0.42)  1(0.93,1.07)  0.36(0.31,0.42) 

-0.13(-0.19,-0.08)  -0.29(-0.34, -0.24)  0.36(0.31,0.42)  1(0.93,1.08) 

By  50  time  units,  the  average  correlation  matrix  stabilizes  and  does  not  change  as  t 

increases;  although,  the  confidence  interval  widths  can  improve  slightly  at  higher  t  values, 
see  Figure  3.3.  The  correlation  matrix  for  N(20)  approaches  the  correlation  matrix  for 
A^(50).  Meanwhile,  the  correlation  matrix  for  A^(l)  provides  a  poor  estimate  with  respect 
to  the  correlation  matrices  at  larger  t  values.  The  underlying  reason  that  the  correlation 
matrix  requires  time  to  stabilize  is  due  to  the  number  of  transitions  that  occur  in  a  given 


(3.14) 
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time.  The  particular  model  for  the  sample  problem  only  has  0.37  expected  transitions  at 
t  =  I,  7.44  expected  transitions  at  t  =  20,  and  18.67  expected  transitions  at  t  =  50.  Based 
on  the  testing  in  this  example,  t  should  be  set  so  that  the  number  of  expected  transitions  is 
at  least  greater  than  twice  the  number  of  model  states. 


State-to-State  Transition 


Figure  3.3:  State-to-State  95%  confidence  interval  widths  for  model  correlation 


Using  Algorithm  1,  the  covariance  of  a  Markov  renewal  process  can  be  reasonably 
estimated.  Assuming  the  data  sample  dictates  the  maximum  value  of  t,  the  modeler  is  only 
required  to  specify  the  number  of  simulation  iterations,  N.  Using  the  estimate  of  the  model 
covariance,  the  goodness  of  fit  of  the  theorized  semi-Markov  process  with  respect  to  the 
observed  data  sample  can  be  determined. 
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IV.  Goodness  of  Fit  Metrics  in  Semi-Markov  Processes 


4.1  The  Goodness  of  Fit  Metric 

The  basic  premise  behind  a  goodness  of  fit  metric  is  to  compare  observed  data  points 
with  their  expectations.  When  the  two  values  are  similar,  the  data  is  said  to  have  a  close  fit 
with  the  model.  In  the  case  of  a  semi-Markov  process,  the  observed  process  value  at  a  time, 
t,  is  the  state  of  the  system  at  t,  while  the  expected  process  value  is  the  average  state  of  the 
system  at  t.  For  a  fixed  time,  t,  a  semi-Markov  process  only  has  one  comparison  available, 
namely  the  actual  occupied  state  versus  the  expected  occupied  state.  This  lone  comparison 
does  not  provide  enough  information  to  determine  the  level  of  similarity  between  the 
expected  and  observed  values,  especially  when  the  expected  occupied  state  is  constant 
once  the  model  reaches  steady-state. 

Consider  a  four  state  semi-Markov  process  that,  on  average,  spends  40%  of  its  time  in 
State  1,  10%  in  State  2,  10%  in  State  3,  and  40%  in  State  4.  The  long  run  expected  value 
for  the  semi-Markov  process  would  be  2.5  which  is  not  an  actual  state  and  is  the  furthest 
possible  value  from  the  two  most  commonly  occuring  states  in  the  system.  In  this  scenario, 
80%  of  the  time  the  absolute  value  of  the  difference  between  expected  and  observed  values 
is  1.5  and  20%  of  the  time  the  absolute  value  is  0.5.  Simply  reordering  the  states  by 
swapping  State  1  with  State  2  and  State  3  with  State  4  would  maintain  the  expectation  of 
2.5,  but  the  absolute  value  of  the  difference  between  expected  and  observed  values  would 
become  0.5  for  80%  of  the  time  and  1.5  for  the  remaining  20%  of  the  time.  This  arbitrary 
change  would  result  in  different  model  fit  estimates  for  two  equivalent  models. 

On  the  other  hand,  as  discussed  in  Section  2.2,  the  semi-Markov  process  possesses  an 
underlying  Markov  renewal  process  that  contains  enough  information  to  determine  the  level 
of  similarity  between  the  observed  data  and  the  model.  At  t,  the  observed  Markov  renewal 
process  is  a  vector  of  the  number  of  visits  to  each  individual  state.  The  corresponding 
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expected  number  of  visits  to  each  state  can  be  calculated  in  the  transform  domain  using  the 
generating  function,  defined  in  Pyke  [15]  as  follows: 

'P,  =  1  -  (1  -  z)(I  -  q)-'q[zl  +  (1  -  z)  A(l  -  q)“M]“'.  (4.1) 


Appendix  A  shows  the  full  derivation  of  the  expected  value  and  the  variance  for  a  Markov 
renewal  process  in  the  transform  domain  based  on  'Pj,.  The  time  domain  expected  value  and 
variance  are  calculated  via  a  Laplace  transform  inversion  using  Equations  3.4  and  3.11. 

Once  a  proposed  semi-Markov  process  is  fit  to  observed  data,  a  sample  covariance 
matrix  may  be  calculated  from  the  observed  and  expected  values  of  the  corresponding 
Markov  renewal  process.  The  covariance  matrix  cannot  be  estimated  using  Equation  3.2 
because  there  is  only  one  N(t)  vector  value  for  any  given  t.  Instead,  the  sample  covariance 
calculation  involves  finding  the  expected  Markov  renewal  process  counts  at  every  transition 
time  during  the  observed  data  run.  These  expected  counts  are  subtracted  from  the  actual 
state  counts  at  their  respective  times.  Assuming  the  semi-Markov  process  contains  n  states 
and  with  N  observed  transitions,  the  result  is  a  series  of  N  vectors,  each  of  length  n.  Then 
the  covariance  of  the  Markov  renewal  process  at  is 

YiO(h)  -  E(4))(0(4)  -  E(4))',  (4.2) 


where  0,(4)  is  the  observed  count  and  E,(4)  is  the  expected  count  of  State  i  at  the  k-th 
transition  time.  Because  the  variances  are  time  dependent,  the  covariance  of  the  Markov 
renewal  process  is  also  time  dependent.  Traditional  goodness  of  fit  testing  based  on 
covariance  matrices  is  not  time  dependent  and  results  in  a  goodness  of  fit  metric,  GoF, 
as  follows: 


l=l  j=l 


(0(i:,v)-E(i;,y))2 

|E(i:,-/)| 


(4.3) 


which  follows  a  Chi-squared  distribution  with  y  =  ^n{n  -l- 1)  degrees  of  freedom.  However, 
the  time  dependency  of  the  Markov  renewal  process  alters  the  GoF  equation  slightly  so 
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that 


GoF(t)  =yy  ^  ^  (4.4) 

|E(i:, 7(0)1 

While  this  seems  like  a  very  minor  change,  as  Section  4.3  will  show,  the  inclusion  of  time 


has  major  ramifications  for  the  GoF(t)  distribution. 


4.2  Baseline  Model 

A  baseline  semi-Markov  process  from  the  sample  problem  of  Section  3.2  will  be  used 
to  illustrate  the  behavior  of  the  GoF{t)  in  Equation  4.4.  Recall  that  the  baseline  semi- 
Markov  process  consisted  of  four  states  that  are  nearly  fully  connected.  A  visualization 
of  the  baseline  semi-Markov  process  appears  in  Figure  4.1.  The  baseline  semi-Markov 
process  was  defined  by  the  transition  probability  matrix 

0  0.75  0  0.25 

0.5  0  0.5  0 

P  =  (4.5) 

0.25  0.45  0  0.3 

0.15  0.35  0.5  0 

and  a  sojourn  distribution  matrix  of 

0  r(2, 2)  0  r(2, 2) 

r(2,i)  0  r(0.5,4)  0 

F  =  ,  (4.6) 

r(0.5,i)  r(0.5,4)  0  r(2,2) 

r(0.5,4)  r(2,2)  r(2,i)  o 

where 

Jnx 

— — (4.7) 

0 

The  result  is  a  baseline  kernel  matrix  of 

0  0.1875te-5  0 

0.5te-‘  0  1.7725r°'^e-?  0 

G(t)  =  .  (4.8) 

0.8862rO'5e-'  1.5952r°'5e-?  0  O.OlSte-'^ 

0.5317r°'5e-?  0.0S15te--2  Q.5te-‘  0 
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Figure  4. 1 :  The  baseline  semi-Markov  proeess  eonneeted  graph 


Using  the  baseline  semi-Markov  proeess,  the  eovarianee  matrix  of  the  related  Markov 
renewal  proeess  at  time  t  is  ealeulated  as 


i:[N(t)]  = 


Var(iVi(0)  Cov(iVi(t)iV2(t))  Cov(Ni(t)N3{t))  Cov(iVi(t)iV4(0) 
Co\iN2(t)Ni{t))  YariNiit))  Cov(iV2(tW3(0)  Cov(iV2(t)iV4(0) 
CowiN^mm  Cowimmit))  Var(iV3(0)  CoviNsHmit)) 
Cov(iV4(0M(0)  Cov(iV4(0^2(0)  Cov(iV4(tW3(0)  Var(iV4(0) 


(4.9) 


Based  on  Equation  3.1,  i:[N{t)]  =  [E[Mi{t)Mj{t)]]  =  [Cov[A^,(t)A^X0]]  =  [E[(M(0  - 
^[Ni(t)])(N j(t)  -  E[A^/t)])]]  whieh  ean  be  simulated.  Thus,  the  resulting  GoF{t)  metrie 
for  a  sample  run  of  the  baseline  model  is 

(OK,(0)  -  E®/!)))" 


GoF(,)  = 

;=i  j=i 


|E(i:,;(0)l 


(4.10) 


Sinee  there  are  ten  elements  along  the  upper  diagonal  of  eaeh  matrix  used  in  the  GoF{t) 
ealeulation,  the  initial  expeetation  is  for  the  GoF{t)  to  follow  a  Chi-squared  distribution 
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with  nine  degrees  of  freedom.  However,  the  inclusion  of  the  time  aspect  actually  results  in 
the  distribution  of  GoFit)  more  closely  following  a  transformed  Chi-squared  distribution. 

4.3  Goodness  of  Fit  Metric  Distribution 

The  distribution  of  GoFit)  is  easily  observed  through  simulation.  The  simulation 
involves  collecting  a  sample  run  of  the  baseline  model  out  to  time  t.  For  each  transition,  k, 
from  time  0  to  t,  the  Markov  renewal  process  count  Nitk)  and  the  E[A/^(4)]  are  calculated, 
where  4  is  the  time  of  the  k-th  transition.  The  observed  covariance  between  States  i  and  j 
for  the  sample  run  becomes 

OiZtjim  =  —  Yj^Niiti)  -  EiNiimNjiti)  -  FiNjiti)).  (4.11) 

^  ^  i=i 

The  expected  covariance  matrix  is  calculated  through  simulation  as  described  in  Section 
3.1.  The  GoFit)  metric  is  found  via  Equation  4.10. 

Eigure  4.2  illustrates  the  distribution  of  GoFit)  for  the  baseline  model  at  t=50  after 
collecting  10,000  samples.  As  a  point  of  reference,  the  Chi-squared  distribution  with  9 
degrees  of  freedom  also  appears  on  the  graph.  While  the  metric  distribution  and  theoretical 
distribution  are  similarly  shaped,  the  metric  definitely  does  not  follow  a  Chi-squared 
distribution  with  9  degrees  of  freedom  due  to  the  shift  along  the  .r-axis.  In  addition, 
the  distribution  of  the  metric  changes  over  time  as  illustrated  in  Eigure  4.3  which  has 
completely  different  jc-axes  in  each  portion  of  the  figure,  attributable  to  the  changing 
distribution. 

Although  the  GoFit)  metric  distribution  differs  for  the  various  runtimes  by  shifting 
along  the  .r-axis,  the  shapes  of  the  distributions  follow  a  general  pattern.  Eigure  4.4 
illustrates  the  upper  and  lower  95%  confidence  interval  cutoffs,  means,  and  medians  for  the 
metric  distributions  at  various  runtimes.  All  four  measures  appear  linear  with  respect  to  the 
model  runtime,  clearly  demonstrating  the  distribution  shift  with  respect  to  time.  Since  the 
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Figure  4.2:  GoF{f)  distribution  for  the  baseline  model  at  time  50  after  10,000  iterations 


GoF{t)  distribution  does  not  follow  a  true  Chi-squared  distribution,  the  testing  in  Chapters 
4-6  utilize  simulation  in  place  of  a  theoretical  distribution  to  evaluate  models. 

4.4  GoF{t)  Metric  Sensitivity 

The  ability  of  the  GoF{t)  metric  to  distinguish  between  changes  in  the  model  is  also 
demonstrated  through  simulation.  This  sensitivity  testing  involves  collecting  data  from  a 
similar,  known  model  and  treating  it  as  if  it  came  from  the  target  model.  To  accomplish 
the  sensitivity  testing,  a  simulation  from  an  altered  version  of  the  baseline  model,  5,  is  run, 
where  at  least  one  value  in  P  or  at  least  one  sojourn  distribution  has  been  changed.  The 
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Figure  4.3:  GoF(t)  distribution  at  various  times  after  10,000  iterations 


altered  model  is  called  the  test  model,  T.  Recall  that  a  semi-Markov  process  is  completely 
defined  by  P  and  the  sojourn  distributions  in  F.  After  running  T  out  to  time  t,  the  Markov 
renewal  process  count  A^t(4)  and  the  E[Ab(4)]  are  calculated  for  the  T  and  B  models  at  the 
time  of  the  k-th  transition  (4).  For  testing  purposes,  the  transitions  in  the  test  model  sample 
run  will  determine  the  individual  4  values.  The  observed  covariance  between  States  i  and 
j  for  the  sample  run  becomes 

1  ^ 

OiZijm  =  ^  -  ECAg./t/)),  (4.12) 

where  Nrjitk)  is  the  Markov  renewal  process  count  for  the  /-th  state  of  the  test  model  at 
time  4  and  A^b,,(4)  is  the  Markov  renewal  process  count  for  the  /-th  state  of  the  baseline 
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Figure  4.4:  GoF(t)  distribution  measures  of  central  tendency  and  dispersion 


model  at  time  4.  The  expected  covariance  matrix  is  still  based  solely  on  the  baseline  model 
from  Equation  3.2  and  GoF(t)  is  found  via  Equation  4.10. 

4.4.1  Simulation. 

In  this  simulation,  the  baseline  model  consists  of  P  defined  in  Equation  4.5  and  the 
sojourn  distributions  defined  in  Equation  4.6.  There  are  four  unique  sojourn  distributions 
(r(2,2),  r(2,l),  r(0.5,4),  and  r(0.5,l))  which  can  be  defined  as  a  vector  of  a  values 
(2, 2, 0.5, 0.5)  and  a  vector  of  p  values  (2, 1,4,1)  respectively.  Eor  testing  purposes,  6 
additional  a  and  /3  vectors  were  created  which  can  be  found  in  Table  4.1.  The  al,  a3,  /32, 
and  /33  vectors  contain  large  changes  to  the  sojourn  distributions  that  the  GoF(t)  metric 
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should  readily  identify  because  the  altered  sojourn  times  will  result  in  much  higher  or 
much  lower  state  visit  counts  than  in  the  baseline  model.  The  remaining  vectors  have 
subtle  sojourn  distribution  changes  that  are  more  difficult  for  the  metric  to  identify  because 
the  state  visit  counts  will  not  be  affected  as  much.  Figure  4.5  illustrates  how  the  sojourn 
distributions  change  with  the  various  parameterizations. 


Table  4.1:  Sojourn  Distribution  Parameters 


a  Vectors 

j3  Vectors 

orl  (Baseline) 

(2,2,0.5,0.5) 

ySl  (Baseline) 

(2,1,4,1) 

Q'2 

(4,4,1,1) 

yS2 

(4,2,8,2) 

Q'3 

(1,1,0.25,0.25) 

yS3 

(1,0.5,2,0.5) 

0-4 

(4,2,0.5,0.5) 

yS4 

(4,1,4,1) 

Q'5 

(2,8,0.5,0.5) 

135 

(2,4,4,1) 

a6 

(2,2,5,0.5) 

/36 

(2,1,40,1) 

al 

(2,2,0.5,50) 

131 

(2,1,4,100) 

Each  combination  of  a-  and  jS-vectors  is  examined  by  creating  a  model  with  those 
specific  parameters.  Each  test  run  simulates  the  model  to  t  =1,000  time  units.  The  test  run 
covariance  matrix  is  created  using  the  Markov  renewal  process  counts  of  the  test  model, 
the  expected  Markov  renewal  counts  of  the  baseline  model,  and  Equation  4.12  for  each 
transition  in  the  test  run.  An  expected  covariance  matrix  is  calculated  using  500,000 
iterations  of  the  baseline  model  and  the  Covariance  Simulation  Algorithm  described  in 
Section  3.1.  The  GoF{t)  metric  for  the  sample  run  is  computed  as  in  Equation  4.10.  This 
process  is  repeated  for  5,000  iterations  to  build  a  distribution  for  the  GoFit)  metric  of  the 
test  model  versus  the  baseline  model. 
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Figure  4.5:  Sojourn  probability  density  funetion  relationships 


To  establish  rejeetion  eriteria,  100,000  iterations  of  the  baseline  model  versus  itself 
were  eondueted  with  the  upper  5%  of  the  test  statistie  values  eonsidered  as  the  rejeetion 
region  for  the  test.  This  rejeetion  region  establishes  a  type  I  error  rate  (ineorreetly  rejeeting 
a  true  model)  of  0.05  for  the  test.  Table  4.2  illustrates  the  type  II  error  rates  (ineorreetly 
failing  to  rejeet  a  false  model)  for  the  test  model  simulations.  In  this  ease,  the  crl-/?!  entry 
is  the  baseline  model  whieh  will  not  have  a  type  II  error,  but  does  have  a  type  I  error  set 
at  0.05.  Of  the  other  48  test  seenarios,  42  seenarios  have  type  II  error  rates  less  than  0.05 
and  only  4  seenarios  have  type  II  error  rates  greater  than  0.25.  These  4  seenarios  are  the 
al  -  /33  model,  the  a3  -  (32  model,  the  al  -  J33  model,  and  the  al  -  [31  model. 
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Table  4.2:  Model  Type  II  Error  Rates 


a  Vectors 

yS  vectors 

al 

a2 

q;3 

aA 

aS 

q;6 

al 

0 

0 

0 

0 

0 

0 

0 

yS2 

0 

0 

0.879 

0 

0 

0 

0 

yS3 

0 

0.973 

0 

0.002 

0.0006 

0.0004 

0.666 

yS4 

0.0002 

0 

0.041 

0 

0 

0 

0 

yS5 

0 

0 

0.023 

0 

0 

0 

0 

(36 

0 

0 

0.068 

0 

0 

0 

0 

yS7 

0.0026 

0 

0.252 

0 

0 

0 

0.062 

Analyzing  the  effects  of  the  a-  and  /3-vectors  explains  the  higher  type  II  error  rates 
in  the  4  non-baseline  models.  For  the  a3  -  /32  model,  the  cr- vector  is  exactly  half  the 
baseline  or- vector  while  the  yS- vector  is  exactly  double  the  baseline  /3- vector.  For  these  F- 
distribution  parameterizations,  the  means  of  the  distributions  equal  a(3  and  the  variances 
equal  a0^.  These  parameterizations  imply  that  the  a\  -/31  and  cr3  - 132  models  possess 
identical  sojourn  distribution  means  but  the  a3-f32  model  has  twice  the  sojourn  distribution 
variances.  The  additional  variance  results  in  the  lower  type  II  error  rate  of  0.879  when 
compared  with  the  a2  - [33  model’s  type  II  error  rate  of  0.973.  The  Qr2  -/33  model  has  half 
the  sojourn  distribution  variances  compared  with  the  baseline  model  while  maintaining  the 
same  means.  In  this  case,  the  test  model  performs  “better”  than  the  true  model  in  that  fewer 
model  runs  will  be  rejected  even  though  the  model  is  wrong. 

The  al  -  /33  and  a3  -  jS7  models  have  type  II  error  rates  greater  than  0.1  due  to  the 
balance  between  the  new  parameterizations.  The  difference  between  these  two  models 
and  the  baseline  model  is  the  r(0.5,l)  distribution  which  becomes  a  r(50,0.5)  distribution 
with  all  other  yS  values  cut  in  half  in  the  al  -  /33  model  and  a  r(0.25,100)  with  all  other 
a  values  cut  in  half  in  the  a3  -  yS7  model.  This  results  in  a  larger  mean  and  variance  for 
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one  distribution  (the  r(50,0.5)  in  the  al  -  /33  seenario  and  the  r(0.25,100)  in  the  a3  -/?7 
seenario)  while  the  other  three  distributions  from  F  in  eaeh  respeetive  seenario  are  redueed. 
The  typieal  run  with  either  of  these  models  has  relatively  few  transitions  due  to  the  large 
mean  of  the  assoeiated  T-distribution.  However,  the  redueed  means  and  varianees  in  the 
other  distributions  eounterbalanee  the  more  misspeeified  distribution. 

Figure  4.6  shows  the  rejeetion  probabilities  for  eaeh  test  model  versus  the  baseline 
model  over  the  time  length  of  the  test  run.  Most  of  the  seenarios  stabilize  with  rejeetion 
rates  above  0.9  by  200  to  500  time  units.  This  equals  a  type  II  error  of  less  than  0.1.  The 
baseline  model  plus  the  five  test  models  diseussed  above  eorrespond  to  the  five  lines  that 
remain  below  0.9  all  the  way  out  to  1,000  time  units.  Naturally,  a  wide  range  of  models 
would  perform  similarly  to  the  baseline  model  as  illustrated  by  these  test  models. 
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Figure  4.6:  Rejeetion  pereentages  for  various  sojourn  distributions  over  time 
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In  addition  to  the  various  sojourn  distribution  parameterizations,  testing  included  a 
variety  of  probability  transition  matrices.  Appendix  B  contains  the  rejection  percentage 
matrices  for  seven  additional  P  matrices.  The  end  result  is  that  this  test  is  sensitive  against 
even  slight  changes  to  the  probability  transition  matrix.  Even  the  a\  - parameterizations 
have  a  much  higher  rejection  rate  when  the  various  P  matrices  are  used.  In  these  cases, 
the  fact  that  the  sojourn  time  distributions  are  identical  is  offset  by  the  differences  in 
the  transition  probabilities.  These  differences  result  in  highly  altered  sample  paths  that 
are  unlikely  to  be  produced  by  the  baseline  model.  In  most  of  the  cases,  the  a?  -  yS? 
parameterization  has  a  slightly  lower  rejection  rate  in  the  80-90%  range  because  the 
misspecified  r(50, 100)  sojourn  distribution  reduces  the  Markov  renewal  process  counts 
to  the  point  that  the  altered  P  is  masked. 

The  reason  that  perturbations  of  P  have  such  a  profound  effect,  in  most  cases,  on  the 
ability  to  distinguish  between  models  is  due  to  the  portions  of  the  covariance  matrix  that 
are  influenced  by  P.  F  effects  the  time  between  transitions.  Suppose  the  distributions  in  F 
were  shifted  along  the  .r-axis  so  that  their  means  were  doubled.  The  impact  of  this  change 
would  be  to  halve  the  Markov  renewal  process  counts  and  expectations  for  a  given  time. 
Ultimately,  the  model  would  behave  similarly  to  the  original  model,  it  would  just  take 
longer  to  transition.  P  effects  the  actual  behavior  of  the  model.  For  the  baseline  model, 
doubling  the  probability  of  transitioning  from  State  1  to  State  4  (from  0.25  to  0.5)  would 
also  reduce  the  probability  of  transitioning  from  State  1  to  State  2  (from  0.75  to  0.5).  Not 
only  does  this  change  alter  the  Markov  renewal  process  counts  and  expectations,  but  it  also 
impacts  the  correlation  matrix  for  the  model  by  reducing  the  correlation  among  States  1  and 
2  and  increasing  the  correlation  between  States  1  and  4.  For  a  given  t,  the  new  correlation 
matrix  results  in  a  new  covariance  matrix  as  well  since  left  hand  side  of  Equation  3.13  has 
changed. 
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The  GoF(t)  metric  developed  in  this  chapter  can  detect  differences  between  the 
observed  data  and  expected  data  for  an  underlying,  hypothesized  semi-Markov  process. 
The  GoF{t)  metric  detected  model  differences  in  92%  of  the  models  tested  in  this  chapter. 
A  natural  extension  of  this  goodness  of  fit  testing  is  to  determine  whether  or  not  a  potential 
covariate  has  an  impact  on  the  hypothesized  model.  The  next  chapter  examines  covariate 
detection  which  impacts  data  collection  requirements  and  overall  model  definition. 
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V.  Testing  a  Semi-Markov  Process  for  the  Presence  of  Covariates 


The  GoF{t)  metric  developed  in  Chapter  4  can  be  used  to  determine  whether  or  not 
a  significant  covariate  is  present  within  a  model.  The  process  involves  calculating  two 
GoF(t)  metric  values  for  a  single  observed  data  sample.  The  first  GoF(t)  value  comes 
from  comparing  the  data  to  a  model  that  contains  the  covariate  and  determining  where 
this  GoF{t)  value  falls  in  the  distribution  of  GoFit)  values  for  this  covariate  model.  The 
second  GoF(t)  value  is  created  by  comparing  the  data  to  a  model  that  does  not  contain  the 
covariate  and  determining  where  this  second  GoF(t)  value  falls  in  the  GoF{t)  distribution  of 
the  non-covariate  model.  The  major  difference  between  the  models  will  be  the  state  space 
representations  which  will  alter  the  expectations  and  covariance  matrices  used  to  create  the 
GoF(t)  metric.  The  changes  in  the  state  space  representation  require  the  covariates  to  be 
categorical  rather  than  continuous,  as  a  continuous  covariate  would  result  in  an  uncountably 
infinite  state  space. 

Before  demonstrating  the  test  for  covariates,  consider  the  baseline  model  to  add 
context.  Suppose  the  four  state  model  in  Figure  4.1  tracks  an  individual  as  he  or  she 
progresses  through  the  various  stages  of  a  disease  where  State  1  =  Healthy,  State  2  = 
Diseased,  State  3  =  Recovering,  State  4  =  Disease  No  Longer  Present.  In  this  scenario,  the 
Healthy  State  only  occurs  prior  to  an  individual  being  diagnosed  with  the  disease  for  the 
first  time.  Any  transitions  back  to  the  Healthy  State  from  the  other  three  states  represent 
the  death  of  the  current  subject  and  the  admission  of  a  new  subject  into  the  system. 

Within  this  scenario,  the  subject’s  sex  could  be  a  covariate  in  the  model  as  disease 
progression  may  substantially  vary  between  men  and  women.  The  eight  state  model  in 
Figure  5.1  represents  a  state  space  expansion  to  account  for  a  subject’s  sex.  In  the  covariate 
model,  let  State  1  =  Healthy  Female,  State  2  =  Diseased  Female,  State  3  =  Recovering 
Female,  State  4  =  Disease  No  Longer  Present  Female,  State  5  =  Healthy  Male,  State  6  = 
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Diseased  Male,  State  7  =  Reeovering  Male,  State  8  =  Disease  No  Longer  Present  Male. 
The  dotted  arrows  represent  transitions  when  a  subjeet  dies  and  the  next  subjeet  is  of  the 
opposite  sex.  Effeetively,  there  are  two  separate  four-state  models  that  are  eonneeted  at 
speeifie  transition  points. 


Figure  5.1:  The  eovariate  baseline  semi-Markov  proeess  eonneeted  graph 


Ultimately,  the  state  spaee  ehanges  required  for  the  eovariate  model  result  in  ehanges 
to  the  transition  probability  matrix,  the  sojourn  distribution  matrix,  and  the  kernel  matrix. 
Let  n  represent  the  probability  that  a  given  subjeet  is  eovariate  level  1  and  \-n  represent 
the  probability  that  a  subjeet  is  eovariate  level  2.  The  starting  eovariate  model  transition 
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probability  matrix  for  each  state,  5,-,  in  the  covariate  baseline,  CB,  model  is 


PcB  - 


Si 

Si 

Ss 

S4 

Ss 

Se 

S7 

Sg 

Si 

0 

0.75 

0 

0.25 

0 

0 

0 

0 

Si 

0.5n 

0 

0.5 

0 

0.5(1  -  n) 

0 

0 

0 

Si 

0.25n 

0.45 

0 

0.3 

0.25(1 -tt) 

0 

0 

0 

S4 

0.15;r 

0.35 

0.5 

0 

0.15(1 -n) 

0 

0 

0 

Ss 

0 

0 

0 

0 

0 

0.75 

0 

0.25 

S6 

0.5n 

0 

0 

0 

0.5(1  -  n) 

0 

0.5 

0 

Si 

0.25n 

0 

0 

0 

0.25(1  -;r) 

0.45 

0 

0.3 

Sg 

0.15;r 

0 

0 

0 

0.15(1  -;r) 

0.35 

0.5 

0 

(5.1) 


(5.2) 


while  the  starting  covariate  model  sojourn  distribution  matrix  from  the  example  for  the  CB 
model  is 

Si  Si  Si  S4  Ss  Se  Si  Sg 

Si  0  Fi  0  Fi  0  0  0  0 

Si  F2  0  F3  0  F2  0  0  0 

Ss  F4  F3  0  Fi  F4  0  0  0 

54  F3  Fi  F2  0  F3  0  0  0 

Fcb  = 

55  0  0  0  0  0  Fi  0 

Sf,  F2  0  0  0  F2  0  F3  0 

Si  F4  0  0  0  F4  F3  0  Fi 

Sg  F3  0  0  0  F3  Fi  F2  0 

where  Fi  is  a  r(2,2),  F2  is  a  r(2,l),  F3  is  a  r(0.5,4),  and  F4  is  a  r(0.5,l).  The  CB  subscript 
for  matrices  is  used  to  denote  the  covariate  baseline  model  matrices  as  opposed  to  the 
original  matrices  used  in  Chapter  4  and  the  5,  notation  tracks  the  states  of  the  covariate 
model.  The  covariate  baseline  kernel  matrix  becomes 

A  B 
C  D 


Gcb  - 


(5.3) 
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where 


0  0.1875;ce-2  0  0.0625xe-i 

O.Snxe-^  0  1.7725;c-0'5e-?  0 

0.8862;r;c-0'5e-^  1.5952;c-°'5e-i  0  Q.Q15xe-i 

Q.53\lnx-^-^e~-^  Qmi5xe-i  Q.Sxe-^  0 

0  0  0  0 

0.5(1  -  7r);ce-^  0  0  0 

B  = 

0.8862(1  -7r);c-°5e-^  0  0  0 

0.5317(1 -7r);c-‘’V?  0  0  0 

0  0  0  0 

0.57r.re“^  0  0  0 

C  = 

0.88627r;c-o5e-^  0  0  0 

Q.53\lnx-^-^e-^^  0  0  0 


0 

0.m5xe--2 

0 

0.0625a:e-5 

0.5(1  -  n)xe~^ 

0 

1.7725a:-°-5e-? 

0 

0.8862(1 

1.5952a:-°'5e-? 

0 

0.015x6-^2 

0.5317(1 

0.0875a:e-i 

0.5xe~’‘ 

0 

A  single  observed  data  set  can  be  tested  against  either  a  four-state  model  or  an  eight- 
state  model  depending  on  whether  or  not  the  covariate  is  taken  into  account.  Consider  a  test 
run  where  four  subjects  are  female,  male,  male,  and  female  respectively.  A  potential  sample 
path  given  these  four  subjects  appears  in  Table  5.1.  If  the  covariate  was  insignificant,  as  is 
the  case  with  the  covariate  baseline  model  depicted  above,  data  collected  from  the  eight- 
state  model  will  pass  the  goodness  of  fit  test  for  the  four-state  model  once  the  data  is 
collapsed  to  four  states.  This  will  be  demonstrated  in  Section  5.1. 
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Table  5.1:  Example  Markov  Renewal  Process  Sample  Paths 


Four-state 

Eight- State 

Markov  Renewal  Process 

Markov  Renewal  Process 

Sample  Path 

Sample  Path 

Subject 

Without  Covariate 

With  Covariate 

1 

1 

1 

2 

2 

3 

3 

1 

5 

2 

4 

8 

2 

6 

1 

5 

3 

2 

6 

3 

7 

1 

1 

4 

4 

4 

3 

3 

A  significant  covariate  would  necessitate  a  change  to  at  least  one  of  the  values  or 
distributions  in  the  Pcb  or  Fcb  matrices  which  also  alters  the  Gcb  matrix.  For  testing 
purposes,  assume  that  the  four-state  baseline  model  does  not  change.  This  assumption 
forces  the  covariate  model  changes  to  be  balanced  within  the  model.  For  instance,  starting 
with  Pcb,  if  the  probability  of  going  from  State  5  to  State  8  is  changed  to  0.5,  then  at  least 
one  of  the  other  transition  probabilities  in  the  State  5  row  must  change  since  rows  must  sum 
to  1.  Assume  that  the  State  5  to  State  6  transition  probability  also  becomes  0.5  to  satisfy 
this  constraint.  In  addition,  the  State  1  to  State  2  and  the  State  1  to  State  4  transitions  would 
change  based  on  n  so  that  their  average  transition  probabilities  when  the  covariate  is  not 
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considered  are  0.75  and  0.25  respectively.  The  new  State  1  to  State  2  and  State  1  to  State  4 
transition  probabilities  for  the  eight-state  model  can  be  calculated  by  solving  the  following 
series  of  equations: 


0.75  =  pnn  -l-  P5(,{1  -  n)  =  pn  *n  +  0.5  *  (1  -  .tt) 

0.25  =  Pun  -I-  j!?58(l  -  n)  =  pu  *  n  +  0.5  *  (1  -  .tt) 

1  =  7*12  +  7*14 
0  <  7’12,7’14,7’56,7’58  <  1 

Generally,  a  covariate  probability  value  n  will  either  be  known  or  assumed,  as  such,  there 
remain  only  three  equations  and  two  unknowns.  This  over- specified  system  of  equations 
limits  the  amount  by  which  any  of  the  individual  probabilities  in  the  Pcb  matrix  may 
change.  In  this  case,  assuming  n  =  0.5,  then  pn  =  1  and  pu  =  0.  If  one  tries  to  make  the 
State  5  to  State  8  transition  probability  0.6  then 

0.75  =  Pu  *  0.5  -I-  0.4  *  0.5 

0.25  =  Pu  *  0.5  -I-  0.6  *  0.5 

1  =  7*12  +  7*14 

implies  pu  =  1.1  and  pu  =  -0.1.  The  probabilities  greater  than  1.0  and  less  than  0.0 
violate  the  laws  of  probability  which  means  this  particular  change  to  the  State  5-State  8 
transition  cannot  occur  in  this  model. 

Similarly,  a  change  to  one  of  the  sojourn  distributions  forces  at  least  one  other  sojourn 
distribution  to  change  if  the  covariate  model  is  going  to  maintain  balance  with  the  smaller 
baseline  model.  In  this  case,  a  balance  exists  between  States  1  and  5  in  the  covariate  model 
versus  State  1  in  the  baseline  model.  Additionally,  States  2  and  6,  States  3  and  7,  and 
States  4  and  8  in  the  covariate  model  are  balanced  with  States  2,  3,  and  4,  respectively, 
in  the  baseline  model.  This  balance  is  based  on  the  assumption  that  the  original  four- 
state  model  provides  an  adequate  fit  of  the  data  and  that  the  larger  covariate  model  is  now 
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under  test.  If  the  covariate  model  distributions  were  not  balanced  with  their  baseline  model 


counterparts,  then  the  covariate  model  would  reduce  to  a  four  state  model  that  is  not  the 
baseline  model.  This  concept  of  balance  is  a  function  of  artificial  semi-Markov  processes 
used  in  this  testing.  In  normal  implementation,  the  data  would  dictate  the  model  parameters 
for  both  the  co variate  model  and  the  reduced  model. 

To  demonstrate  the  sojourn  distribution  balance,  if  the  sojourn  time  for  transitioning 
from  State  5  to  State  8  followed  a  r(3,2)  distribution,  then  the  transition  distribution  from 
State  1  to  State  4  in  the  covariate  model  would  also  have  to  change  to  ensure  the  average 
sojourn  distribution  between  the  two  transitions  remains  a  r(2,2)  distribution  which  is  the 
sojourn  distribution  between  States  1  and  4  in  the  baseline  model.  While  there  are  a  variety 
of  ways  to  balance  the  distributions,  for  simplicity,  the  distributional  mean  is  being  kept 
consistent.  In  the  case  of  a  r(2,2)  distribution  changing  to  a  r(3,2)  distribution,  the  State 
1  to  State  4  sojourn  distribution  would  become  a  r(l,2)  distribution  assuming  n  =  0.5. 
Using  this  pair  of  distributions,  the  average  of  all  State  1  to  State  4  and  all  State  5  to 
State  8  transitions  is  4  time  units,  matching  the  original  mean  of  the  r(2,2)  distribution. 
This  method  does  not  maintain  the  distributional  variance  as  the  combination  of  the  r(l,2) 
distribution  and  the  r(3,2)  distribution  has  a  variance  of  12  while  the  original  variance  of 
the  r(2,2)  distribution  is  8.  In  reality,  the  sample  data  itself  would  be  used  to  determine  the 
P  and  F  matrices  for  the  model  with  covariates  and  the  model  without  covariates. 

5.1  Covariate  Testing 

Once  the  two  models  (baseline  and  covariate)  are  selected,  testing  proceeds  by 
calculating  the  GoF(t)  metric  of  the  sample  data  for  both  models  using  the  methods 
outlined  in  Chapter  4.  The  GoFit)  values  indicate  whether  or  not  either  of  the  models 
can  be  rejected  as  poor  fits  for  the  data.  In  the  case  where  both  models  fail  to  reject,  a 
preference  must  be  made  with  respect  to  either  the  larger  covariate  model  or  the  model 
without  covariates.  This  decision  is  ultimately  up  to  the  modeler  who  must  balance  data 
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collection  requirements  and  computational  complexity  with  the  added  model  fidelity  that 
the  covariates  provide. 

In  this  testing,  the  baseline  four-state  model  from  Chapter  4  is  the  model  without 
covariates  and  a  suite  of  models  with  covariates  was  created  to  illustrate  the  ability  of 
the  GoF(t)  metric  test  to  differentiate  between  the  models.  The  first  round  of  testing 
demonstrates  the  equivalence  between  the  four-state  baseline  model  from  Chapter  4  and 
eight-state  covariate  baseline  model  defined  by  Pcb  and  Fcb-  The  testing  involves  drawing 
A^=5,000  iterations  from  the  covariate  baseline  model  using  n  =  0.5.  The  eight-state  sample 
paths  are  reduced  down  to  four  states  by  making  all  State  5  instances  in  the  sample  path 
State  I’s,  all  State  6’s  instances  State  2’s,  all  State  7’s  instances  State  3’s,  and  all  State 
8’s  instances  State  4’s.  This  transition  effectively  ignores  the  potential  presence  of  the 
covariate.  Ignoring  the  covariate,  4.58%  of  the  runs  are  rejected  based  on  their  GoFit) 
metric  values  which  is  in  line  with  the  expected  5%  that  should  be  incorrectly  rejected.  For 
all  of  the  covariate  testing,  a  type  I  error  rate  of  5%  is  used. 

Ideally,  any  covariate  model  that  is  properly  balanced  with  respect  to  a  four-state 
model  will  have  a  rejection  rate  around  5%  when  it  is  reduced  to  the  four-state  model 
without  covariates.  The  difference  in  models  can  be  detected  by  comparing  the  larger 
covariate  models  to  one  another.  The  first  step  is  a  demonstration  that  the  baseline  model 
and  covariate  baseline  model  are  equivalent  to  one  another  in  the  larger  eight-state  domain. 
This  round  of  testing  draws  /7=5,000  iterations  from  the  four-state  baseline  model  with  one 
slight  difference.  Each  time  the  simulation  returns  to  State  1,  signifying  a  new  subject  enters 
the  system,  the  simulation  determines  whether  or  not  the  new  subject  is  male  or  female 
according  to  covariate  level  probability,  n.  After  the  simulation  concludes,  the  Markov 
renewal  process  is  expanded  to  eight  states  by  making  States  1,  2,  3,  and  4  into  States  5,  6, 
7,  and  8  if  a  particular  subject  is  male.  Note  that  for  an  assumed  tt  =  0.5,  Pij  =  p(,+4)(y+4)  and 
Fij  =  F(,+4)(j+4)  for  1  <  /,  j  <  4.  GoF(t)  metric  values  are  then  calculated  for  the  runs  by 
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comparing  them  directly  with  the  eight-state  covariate  baseline  model.  In  this  case,  5.04% 
of  the  samples  are  rejected  which  is  in  line  with  the  expected  5%  that  should  be  incorrectly 
rejected. 

From  this  point  forward,  any  potential  covariate  model  is  tested  directly  against  the 
covariate  baseline  model.  The  underlying  analysis  is  that  if  the  data  is  not  significantly 
different  from  the  covariate  baseline  model,  the  proposed  data  model  can  be  reduced  down 
to  the  four-state  baseline  model  without  losing  any  fidelity  at  the  modeler’s  discretion. 
Practically,  this  means  it  is  not  necessary  to  track  the  presence  of  the  covariate  because  it 
does  not  actually  add  any  additional  information  about  the  process. 

Testing  focused  on  two  separate  probability  transition  matrices.  Pci  and  Pc2,  and  four 
different  sets  of  Gamma-distribution  parameters,  a2  0:3  -1^3,  ctA  -  A,  and  -/Ss.  The 
first  probability  transition  matrix. 


Si  S2  S3  S4  Ss  Se  Si  Ss 

s,  0100  0  000 

52  0.5n  0  0.5  0  0.5(1  -;r)  0  0  0 

53  0.25;r  0.45  0  0.3  0.25(1  -;r)  0  0  0 

54  0.15;r  0.35  0.5  0  0.15(1  -;r)  0  0  0 

Pci  =  ,  (5.4) 

55  0  0  0  0  0  0.5  0  0.5 

56  O.Stt  0  0  0  0.5(1  -;r)  0  0.5  0 

5,  0.25;r  0  0  0  0.25(1  -;r)  0.45  0  03 

58  0.15;r  0  0  0  0.15(1  -;r)  0.35  05  0 
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represents  a  drastic  change  in  comparison  with  the  covariate  baseline  model  for  the 
probability  transitions  of  States  1  and  5.  The  second  probability  transition  matrix, 

Si  52  53  54  55  56  5,  58 

5.  0  0.875  0  0.125  0  0  0  0 

52  0.5;r  0  0.5  0  0.5(1  -;r)  0  0  0 

53  0.25;r  0.45  0  0.3  0.25(1  -tt)  0  0  0 

54  0.15;r  0.35  0.5  0  0.15(1  -tt)  0  0  0 

Pc2  =  ,  (5.5) 

55  0  0  0  0  0  0.625  0  0.375 

56  0.5;r  0  0  0  0.5(1  -;r)  0  0.5  0 

57  0.25;r  0  0  0  0.25(1  -tt)  0.45  0  03 

58  0.15;r  0  0  0  0.15(1  -tt)  0.35  05  0 

represents  a  smaller  change  for  the  probability  transitions  of  States  1  and  5.  For  the  F- 
distribution  parameters,  let  the  upper  and  lower  halves  of  the  Fcb  matrix  each  have  four 
separate  distributions  (F(2,2),  F(2,l),  F(0.5,4),  F(0.5,l)).  Following  the  notation  used  in 
Section  4.4.1,  define  a  vector  of  a  and  jS  values  for  all  eight  states  as  (2,2,0.5,0.5,2,2,0.5,0.5) 
and  (2, 1,4, 1,2, 1,4,1)  respectively.  Table  5.2  shows  the  original  baseline  a  and yS  vectors  and 
the  four  additional  sets  of  vectors  tested.  Of  the  four  test  vectors,  the  first  three  only  change 
a  single  pair  of  a  or  jS  values,  while  the  fourth  set  represents  a  more  drastic  change  across 
all  distributions.  The  distributions  were  selected  to  balance  the  baseline  distributions  with 
respect  to  the  combined  average  of  each  pair  of  sojourn  distributions  when  compared  with 
the  same  pair  in  the  covariate  baseline  model.  Basically,  the  average  among  all  State  1  to 
State  2  transitions  and  all  State  5  to  State  6  transitions  is  the  same  whether  the  covariate 
baseline  vectors  or  any  of  the  four  additional  vector  pairs  are  used.  The  need  to  maintain 
this  balance  eliminates  the  need  to  run  the  combinatorial  analysis  with  every  P  -  a  -  /? 
permutation  as  seen  in  Chapter  4. 
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Table  5.2:  Covariate  Sojourn  Distribution  Parameters 


a  Veetors 

jS  Veetors 

al  (Baseline) 

(2, 2, 0.5, 0.5, 2, 2, 0.5, 0.5) 

/31  (Baseline) 

(2,1,4,1,2,1,4,1) 

a'2 

(3,2,0.5,0.5,12,0.5,0.5) 

yS2 

(2,1,4,1,2,1,4,1) 

Q'3 

(2,2,0.5,0.5,2,2,0.5,0.5) 

yS3 

(2,1,1, 1,2,1,7,1) 

cr4 

(60,2,0.5,0.5,0.1,2,0.5,0.5) 

j84 

(0.1,1,4,1,20,1,4,1) 

cr5 

(3,1,1,0.25,3,3.5,0.5,0.75) 

j85 

(1,0.5,3,0.25,1.67,1,2,1.25) 

Comparing  the  VciaiPi  model  to  the  baseline  eovariate  model  (PcBCtiPi),  the  type 
II  error  rate  is  0.499  while  the  VcifXiPi  model  has  a  type  II  error  rate  of  0.891  when 
eompared  with  the  baseline  eovariate  model.  The  disparity  in  type  II  error  rates  is  expeeted 
beeause  the  PciOfiySi  model  represents  a  larger  ehange  to  P  than  the  Vc2CtiPi  model.  For  the 
^CBCtiPi,  PcbQ'sAj  ^cbc^aPa,  and  PcBcrj^Ss  models,  the  type  II  error  rates  when  eompared 
with  the  baseline  eovariate  model  are  0.934,  0.933,  0.392,  and  0.948,  respeetively.  Figure 
5.2  shows  the  various  distributions  used  in  the  eovariate  baseline  model  and  the  VcBCtiPi 
and  PceCTsySa  models.  In  both  graphs,  the  solid  line  represents  the  distribution  used  in  the 
eovariate  baseline  model  while  the  two  dashed  lines  represent  the  balaneed  distributions 
from  the  test  models.  The  large  degree  of  overlap  among  the  various  distributions  results 
makes  it  impossible  for  the  GoF{t)  metrie  to  diseem  between  these  models.  On  the  other 
hand,  the  distributions  used  in  the  VcbC(aPa  model  differ  signifieantly  from  the  eovariate 
baseline  model  as  illustrated  in  Figure  5.3.  These  differenees  are  pieked  up  in  over  60%  of 
the  models  tested. 

While  it  may  appear  that  eovariate  testing  laeks  power  to  differentiate  between  models 
based  on  some  of  the  seenarios,  the  reality  is  that  it  ean  deteet  signifieant  differenees  sueh 
those  shown  in  Figure  5.3.  The  seenarios  used  in  this  ehapter  were  speeifieally  erafted  to 
ensure  the  average  aeross  the  eovariate  would  eollapse  down  to  the  baseline  model.  The 
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aipi  vs.  a2p2 


a1  pi  vs.  a3p3 


y  y 

Figure  5.2:  Distributions  from  the  baseline  eovariate,  02/^2,  and  models 


aipi  vs.  a4p4 


Figure  5.3:  Distributions  from  the  baseline  covariate  and  a4/34  models 


end  result  in  several  scenarios  was  a  small  enough  difference  that  the  covariate  proved 
insignificant.  For  an  actual  data  sample,  the  baseline  and  covariate  models  would  be 
dictated  by  the  sample  and  the  underlying  theory  which  results  in  more  realistic  differences 
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between  models.  To  this  end,  the  next  chapter  examines  a  real  world  application  of  the 
GoFit)  metric  with  an  emphasis  on  covariate  detection. 
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VI.  BMI  Model  Analysis 


According  to  the  Centers  of  Disease  Control  and  Prevention  (CDC)  website  [24], 
childhood  obesity  affects  17%  of  US  children  between  the  ages  of  2  and  19,  with  20.5% 
of  US  children  between  the  ages  of  the  12  and  19  considered  obese  as  of  2011  based  on 
Body  Mass  Index  (BMI).  Effectively,  one  out  of  every  five  American  teenagers  is  above 
the  95th  percentile  of  their  BMI-for-age  growth  chart.  A  recent  study  using  regression 
techniques  [25]  concluded  that  obesity  trends  began  rising  in  the  1980s  while  a  study 
using  prevalence  rates  [26]  identified  an  increasing  trend  over  a  span  of  merely  four  years. 
Previous  studies  have  drawn  correlations  between  childhood  obesity  and  various  childhood 
and  adult  afflictions  including  heart  disease  [27]  and  diabetes  [28].  Due  to  these  life- 
threatening  potential  consequences,  many  social  initiatives  have  started  in  recent  years 
targeting  poor  dieting  and  lack  of  exercise,  two  of  the  primary  causal  factors  in  childhood 
obesity  according  to  the  CDC. 

Whether  it  is  the  National  Football  League’s  Play  60  exercise  campaign  or  the 
Healthy,  Hunger-Free  Kids  Act  (the  school  lunch  reform  initiative  championed  by  First 
Lady  Michelle  Obama),  analysts  and  policy  makers  require  robust  tools  to  determine 
if  the  various  initiatives  are  impacting  the  childhood  obesity  rates.  Quickly  identifying 
improvements  in  obesity  rates  can  lead  to  increased  focus  and  funding  for  current 
initiatives,  while  identifying  a  lack  of  improvement  can  result  in  renewed  efforts  to  create 
new  initiatives  to  combat  obesity.  However,  the  community  requires  a  tool  that  can 
differentiate  between  improving  obesity  rates,  stagnating  obesity  rates,  escalating  obesity 
rates  and  spurious  fluctuations  in  the  data  or  noise  that  can  lead  to  false  conclusions.  To  this 
end,  the  GoF(t)  metric  outlined  in  Equation  4.3  can  help  determine  the  tempo  of  current 
childhood  BMI  trends.  While  the  aforementioned  studies  focus  on  whether  or  not  the 
obesity  trends  appear  to  be  changing  over  time,  the  GoF(t)  metric  can  help  determine 
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whether  the  ehanges  are  possibly  random  noise,  merely  a  eonsisteney  in  assoeiation,  versus 
a  true  population  shift. 

Data  from  the  Pels  Longitudinal  Study,  whieh  eontains  the  age  and  BMI  values  for 
ehildren  from  the  1930s  to  the  present,  was  examined  to  explore  the  population  trends  and 
transitions  to  obesity  in  ehildren  over  time.  This  data  eonsisted  of  repeated  measurements 
of  1 146  ehildren  (583  boys  and  563  girls)  at  6  month  intervals  from  age  2  through  age  19  for 
a  total  of  10,771  boy  observations  and  10,434  girl  observations.  On  average,  the  boys  spent 
11.27  years  in  the  study  while  the  girls  spent  1 1.38  years  in  the  study.  Details  regarding  the 
history  and  methods  of  the  Pels  Longitudinal  Study  ean  be  found  in  Roehe  [29].  Based  on 
the  BMI-for-age  growth  eharts  and  the  traditional  BMI-based  eategories  of  Underweight 
(less  than  the  5th  pereentile).  Normal  (5th  to  85th  pereentile).  Overweight  (85th  to  95th 
pereentile),  and  Obese  (greater  than  the  95th  pereentile),  eaeh  ehild  is  traeked  as  he  or  she 
transitions  from  one  eategory  to  the  next.  Pigure  6.1  shows  a  simple  state  model  of  the  four 
BMI-based  eategories.  The  solid  arrows  represent  moving  among  adjoining  states  while 
the  dashed  arrows  indieate  a  jump  of  two  or  more  states. 


Pigure  6.1:  Pour  state  weight  elassifieation  model 


While  skipping  a  state  generally  does  not  oeeur  in  real  life  without  extenuating  faetors 
like  weight  loss  surgery  or  growth  hormones,  the  modeling  proeedure  and  the  six  month 
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(discrete  time)  intervals  in  data  collection  allow  this  phenomenon  to  occur  in  this  model. 
For  this  model,  a  child  enters  the  system  at  the  category  of  their  first  data  measurement. 
The  model  changes  states  when  the  child’s  category  changes,  as  calculated  from  subsequent 
measurements.  After  the  child’s  final  visit,  a  new  child  enters  the  system  at  his  or  her 
current  category  state  and  the  system  proceeds.  If  the  new  child’s  starting  state  is  the 
same  as  the  previous  child’s  ending  state,  the  system  will  continue  accruing  time  in  that 
particular  state  without  transitioning.  In  the  data  set,  this  occurs  227  times  among  boys  and 
197  times  among  girls.  On  the  other  hand,  if  the  new  child  starts  as  Underweight  and  the 
previous  child  ended  up  as  Obese,  the  system  would  undergo  a  transition  from  Obese  to 
Underweight,  skipping  Overweight  and  Normal  weight  in  the  process.  If  a  child  only  had 
one  measurement,  that  child  was  removed  from  further  analysis. 

Based  on  the  observation  time  frames  of  the  children,  two  separate  models  were 
created  from  the  dataset.  The  1930-1960  Model  includes  all  of  the  measurements  prior 
to  1960  and  the  1980-2007  Model  contains  every  measurement  after  1979.  In  the  1930- 


1960  Model, 


P  1930-1960Morfe/ 


uw 

N 

ow 

0 

uw 

0 

0.992 

0.008 

0 

N 

0.338 

0 

0.622 

0.040 

ow 

0.007 

0.850 

0 

0.143 

o 

0.019 

0.264 

0.717 

0 

(6.1) 
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and 


Fl930-1960Mode;  - 


UW 

N 

OW 

0 

UW 

0 

F(0.36,4.82) 

F(0.20,4.00) 

0 

N 

F(0.41, 15.77) 

0 

F(0.43, 19.38) 

F(0.69,7.66) 

OW 

F(0.20,4.00) 

F(0.60, 1.91) 

0 

F(  1.75, 1.02) 

0 

F(0.20,4.00) 

F(0.40,4.98) 

F(0.40, 2.92) 

0 

where  UW  is  Underweight,  N  is  Normal,  OW  is  Overweight,  and  O  is  Obese.  For  the 
1980-2007  Model, 


and 


P  1980-2007Morfe/ 


UW 

N 

OW 

0 

UW 

0 

0.969 

0.031 

0 

N 

0.364 

0 

0.547 

0.089 

OW 

0.016 

0.607 

0 

0.377 

0 

0.032 

0.411 

0.557 

0 

(6.3) 


Model  - 


UW 

N 

OW 

0 

UW 

0 

F(0.76, 2.80) 

F(0.25,5.00) 

0 

N 

F(0.45, 19.90) 

0 

F(0.66, 17.82) 

F(1.51,8.34) 

OW 

F(0.20,4.00) 

F(  1.08, 1.62) 

0 

F(  1.42, 1.44) 

0 

F(0.64,8.58) 

F(  1.09, 3.40) 

F(0.87, 2.34) 

0 

From  a  modeling  standpoint,  the  impaet  of  inereasing  obesity  rates  would  be  most 
notieeable  in  two  main  places,  the  last  column  of  P  and  the  last  row  of  F.  The  final  P 
column  contains  the  probabilities  of  transitioning  into  the  Obese  category  from  one  of  the 
other  categories.  If  obesity  rates  increase,  one  or  more  of  these  transition  probabilities 
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could  increase  which  would  indicate  that  a  higher  portion  of  the  population  reaches  the 
Obese  state.  Meanwhile,  the  bottom  row  of  F  contains  the  distributions  of  the  sojourn 
times  for  the  system  to  leave  the  Obese  state  once  it  arrives  there.  Increases  in  sojourn 
times  for  the  Obese  state  would  indicate  that  the  average  child  that  enters  this  state  spends 
longer  there  and  that  there  is  a  greater  chance  for  the  next  child  to  enter  the  system  in  the 
Obese  category  which  would  not  require  a  state  change. 

Comparing  the  two  models  above,  the  probability  of  reaching  the  Obese  state  is  lowest 
for  the  1930-1960  Model  and  highest  for  the  1980-2007  Model.  From  Table  6.1,  the  1980- 
2007  Model  has  the  largest  average  sojourn  times  for  every  state.  Figure  6.2  shows  the 
difference  in  the  sojourn  time  distribution  going  from  State  2  to  State  4  in  the  two  models. 
The  1930-1960  Model  has  typically  shorter  sojourn  times  for  this  transition  when  compared 
with  the  1980-2007  Model.  Taken  together,  the  1980-2007  Model  implies  that  a  larger 
portion  of  the  child  population  reaches  the  Obese  category  and  on  average  stays  there  for 
a  longer  period  of  time.  The  larger  question  is  whether  the  differences  in  the  models  are 
significant.  The  1930-1960  and  1980-2007  Models  were  tested  against  each  other  after 
t=1000  person  years,  approximately  5  actual  years  of  data  collection.  For  this  study,  a 
person  year  is  the  collective  amount  of  time  the  study  participants  age  during  observation 
which  is  indicative  of  the  number  of  measurements  recorded.  For  instance,  tracking  20 
individuals  for  5  years  would  contribute  100  person  years  to  the  study  which  corresponds 
to  approximately  200  individual  measurements.  The  data  set  contained  1146  children  and 
12,980  person  years.  The  1930-1960  Model  has  a  type  II  error  of  0.441  against  the  1980- 
2007  Model.  This  type  II  error  rate  deceases  to  0.0 1  after  2000  person  years  (roughly 
10  years  of  data  collection).  The  actual  1930-1960  data  sample  is  also  rejected  against 
the  1980-2007  Model.  The  1980-2007  Model  has  a  type  II  error  rate  of  0.076  against  the 
1930-1960  Model  after  1000  person  years  and  the  actual  1980-2007  data  sample  is  rejected 
when  compared  with  the  1930-1960  Model. 
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-  Gamma(0.69,7.66) 

-  Gamma(1.51,8.34) 
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Figure  6.2:  1930-1960  F24  distribution  versus  the  1980-2007  F24  distribution 


Table  6.1:  Average  State  Sojourn  Time  (Years) 


1930-1960  Model 

1980-2007  Model 

Underweight 

1.72 

2.10 

Normal 

7.58 

10.81 

Overweight 

1.23 

1.85 

Obese 

1.38 

2.83 

The  disparity  in  the  type  II  error  rates  is  due  to  the  behavior  of  the  GoF{t)  metrie 
distribution  for  the  two  model  eomparisons.  Speeifieally,  over  40%  of  the  1930-1960 
Model  GoF{t)  metrie  values  still  appear  to  be  from  the  1980-2007  Model.  In  these  eases, 
there  are  enough  visits  to  the  Overweight  and  Obese  States  to  mimie  the  1980-2007  Model 
behavior.  As  t  inereases,  the  1930-1960  Model  GoF{t)  metrie  distribution  separates  itself 
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from  the  1980-2007  Model  GoF{t)  metric  distribution  and  the  type  II  error  rate  decreases 
to  0.01.  The  bottom  line  is  there  exists  strong  evidence  that  the  underlying  BMI  model  has 
changed  significantly  from  the  mid- 1900s  to  the  late- 1900s. 

It  is  worth  noting  that  the  changing  P  and  F  matrices  between  the  1930-1960  Model 
and  the  1980-2007  Model  offer  a  stark  contrast  with  the  matrices  from  the  scenarios  in 
Chapter  5.  While  the  individual  changes  to  P  and  F  are  relatively  small  between  the  two 
BMI  models,  the  cumulative  effect  is  two  models  that  are  relatively  easy  to  distinguish 
between.  In  the  Chapter  5  scenarios,  each  test  model  was  carefully  created  to  ensure  a 
balance  with  respect  to  the  baseline  model,  with  the  added  result  that  competing  models 
were  often  difficult  to  differentiate  between.  This  is  a  case  where  the  real  data  example 
proved  to  be  an  easier  demonstration  of  the  GoFit)  metric  than  the  simulated  scenario,  but 
simplifying  assumptions  regarding  the  modeling  of  transitions  between  states  from  child 
to  child  affect  the  results.  Further  model  refinement  would  be  necessary  to  remove  these 
effects  entirely  for  less  biased  estimates  of  the  transition  probabilities  and  distributions 
between  these  two  cohorts  of  children. 

For  analysts  and  policy-makers  the  utility  of  this  goodness  of  fit  metric  comes  from 
the  ability  to  determine  if  a  data  sample  could  be  the  product  of  a  given  model.  In  the 
case  of  childhood  obesity,  the  metric  can  test  whether  current  obesity  trends  differ  from 
previous  generations.  Suppose  analysts  are  interested  in  the  effectiveness  of  Play  60,  the 
Healthy,  Hunger-Free  Kids  Act,  and  all  of  the  other  recent  programs  targeting  obesity.  After 
collecting  five  years  worth  of  pediatric  visits  for  a  group  of  patients,  the  summary  statistics 
can  tell  whether  the  obesity  rates  are  decreasing.  However,  comparing  this  data  to  the  1980- 
2007  Model  can  tell  whether  there  is  an  actual  change  in  population  behavior  that  warrants 
a  new  model,  or  if  the  current  trends  are  merely  the  product  of  the  stochastic  noise  that  is 
expected  based  on  the  1980-2007  Model.  In  the  former  case,  policy-makers  may  elect  to 
continue  with  the  current  activities,  while  large-scale  program  changes  may  be  warranted 
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in  the  later  case.  In  either  scenario,  policy-makers  have  an  improved  grasp  of  what  the 
current  trends  mean  with  respect  to  the  previous  childhood  obesity  model  behavior. 
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VII.  Conclusion 


Assessing  the  fit  of  a  multi-state  model  with  respeet  to  the  observed  data  used  to  ereate 
the  model  is  an  important  step  in  the  diagnosties  phase  of  modeling.  Previous  diagnosties 
tools  using  prevalenee  eounts  and  Chi-squared  goodness  of  fit  metries  tended  to  possess 
limited  utility  when  eompared  with  the  wide  range  of  models  and  data  eolleetion  teehniques 
used.  The  Goodness  of  Fit  metrie  outlined  in  this  dissertation  has  been  shown  to  overeome 
the  limitations  of  previously  available  diagnostie  tools  by  providing  a  straightforward 
method  of  relating  the  observed  data  direetly  to  a  hypothesized  multi-state  model. 

While  validating  the  goodness  of  fit  for  a  model  is  an  important  diagnostie  step, 
deeiding  among  eompeting  models  is  the  ultimate  goal  of  model  building  with  respeet 
to  providing  a  deeision  maker  with  the  best  possible  information.  To  this  end,  this  work 
illustrates  how  models  that  eontain  eovariate  faetors  ean  be  tested  against  one  another  to 
determine  whieh  model  best  represents  the  true  observed  data.  Proper  model  seleetion  ean 
impaet  model  size,  model  eomplexity,  and  data  eolleetion  efforts  whieh  in  turn  ean  have  a 
profound  impaet  on  the  time  required  to  analyze  a  partieular  model. 

Future  work  may  explore  the  impaet  of  sojourn  distributions  that  laek  a  linear  Laplaee 
transformation  funetion.  In  these  eases,  the  Laplaee  transform  inversion  teehnique  used 
to  ealeulate  expeeted  values  beeomes  a  numerieal  integration  of  a  numerieal  integration 
of  a  funetion  rather  than  the  single  numerieal  integration  of  a  funetion  demonstrated  in 
this  work.  Additionally,  more  work  is  required  in  fitting  the  tails  of  the  GoFiT)  metrie’s 
distribution  with  respeet  to  the  transformed  Chi-squared  distribution.  The  ability  to  test  an 
observed  data  sample  against  a  theoretieal  distribution  rather  than  a  simulated  distribution 
would  greatly  reduee  the  eomputational  workload  in  determining  whether  or  not  a  model 
is  likely  to  produee  an  observation. 
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Appendix  A:  Markov  Renewal  Process  Moments 


The  covariance  calculations  in  Chapter  3  require  the  expected  value  and  variance  of  a 
Markov  renewal  process  at  a  specific  time,  t.  The  following  is  a  derivation  of  both  values 
in  the  transform  domain  based  on  the  generating  function 

¥,  =  !-(!-  z)(I  -  q)-^q{zl  +  (1  -  z)  ^{(I  -  (AT) 

defined  by  Pyke  in  Corollary  5.1  [15].  In  this  case,  z  is  a  random  value  between  0  and  1 
and  q  =  q(t)  =  P  o  ^F(t).  Pyke  uses  'Pj.  to  derive  £(E[N(t)])  in  Theorem  5.2  [15].  Below 
is  a  derivation  of  £(E[N(t)]),  X(E[A(0^]),  and  Var[A(0]- 


£(E[Nm  = 

oz 


z=l 


Consider 


Let  c  = 


d_ 

dz 


X(E[A(0])  = 


-1  -  (1  -  z)(I  -  qy^qVzl  +  (1  -  Z)  T(I  - 
oz 

d{(i--qy^]  ,then 

(I  -  qy^qVzl  +  (1  -  z)c]~^  -  (1  -  z)(I  -  qy^q{-l)* 

[zl  +  (1  -  z)c]“Hl  -  c][zl  +  (1  -  z)c]“^ 

(I  -  qy^q[z\  +  (1  -  z)c]~^  +  (1  -  z)(I  -  qy^q[z\  +  (1  -  z)c]~^* 
[I  -  c][zl  +  (1  -  z)c]“^  and 

(I  -  qy^q[z\  +  (1  -  z)c]~^  +  (1  -  z)(I  -  qy^q[z\  +  (1  -  z)c]“^* 
[I-c][zI  +  (l  -z)c]~^|^^j . 

(I  -  qy\. 


(A.2) 
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~  o  (9^  d 

£(E[N(tf])  =  . 

oz  l  oz  z=l 

d 

=  ^(I  -  gT^gIzI  +  (1  -  Z)c]“^  +  (1  -  z)(I  -  gT^G* 
oz^  oz 

[zl  +  (l  -z)c]“Hl  -  c][zl  +  (1  -z)cY^. 

=  (-1)(I  -  gT^gVzI  +  (1  -  z)c]~^(I  -  c)[zl  +  (1  -  z)c]~^- 
(J  -  qy^ q[zl  +  {I  -z)c]“Hl  -  c][zl  +  (1  -z)c]~^  + 

(1  -z)(I  -  <2')“V(-l)[zI  +  (1  -z)c]“Hl  -  c][zl  +  (1  -z)c]“^* 
[I  -  c][zl  +  (1  -  z)c]“^  +  (1  -  z)(I  -  qT^qVzl  +  (1  -  z)c]“^* 

[I  -  c](-l)[zl  +  (1  -  z)c]“'[I  -  c][zl  +  (1  -  z)c]~^ 

=  -2(1  -  q)~^q[zl  +  (1  -  z)c]~Vl  -  c)[zl  +  (1  -  z)c]~^- 

2(1  -  z)(I  -  qy^qizl  +  (1  -  z)c]~^[I  -  c][zl  +  (1  -  z)c]“^* 

[I  -  c][zl  +  (1  -  z)c]“\  then 
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2(1  -  z)(I  -  qy^qizl  +  (1  -  z)c]~^[I  -  c][zl  +  (1  -  z)c]“^* 
[I-c][zI  +  (l  -z)c]“^|^^j . 
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X(E[^(0'])  =  -2(1  -  q)-^q{l  -  c)  +  (I  -  qT^q. 

Var[^(0]  =  £-\Emf])  -  (£-\E[Nmf  ■ 
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Appendix  B:  Probability  Transition  Matrix  Test  Results 


Using  the  test  eonstruet  from  Chapter  4,  GoF{t)  metrie  tests  utilizing  differing  P- 
matrices  were  also  eondueted.  Eaeh  table  below  eorresponds  to  one  of  seven  P-matriees. 


Pi  = 


0 

0.5 

0.25 

0.15 


0.75 

0 

0.45 

0.35 


0.25 

0 

0.3 

0 


0 

0.5 

0 

0.5 


P2 


0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0.65 

0.10 

0.25 

0.4 

0 

0.4 

0.2 

0.25 

0.45 

0 

0.3 

0.15 

0.35 

0.5 

0 

0 

0.25 

0 

0.75 

0.5 

0 

0.5 

0 

0.45 

0.25 

0 

0.3 

0.5 

0.35 

0.15 

0 

0 

0.15 

0 

0.85 

0.5 

0 

0.5 

0 

0.25 

0.45 

0 

0.3 

0.15 

0.35 

0.5 

0 

(B.l) 


(B.2) 


(B.3) 


(B.4) 


(B.5) 
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0 

0.75 

0 

0.25 

0.05 

0 

0.95 

0 

0.25 

0.45 

0 

0.3 

0.15 

0.35 

0.5 

0 

0 

0.75 

0 

0.25 

0.5 

0 

0.5 

0 

0.75 

0.15 

0 

0.1 

0.15 

0.35 

0.5 

0 

0 

0.75 

0 

0.25 

0.5 

0 

0.5 

0 

0.25 

0.45 

0 

0.3 

0.145 

0.85 

0.005 

0 

(B.6) 


(B.7) 


(B.8) 


Table  B .  1 :  Type  II  Error  For  P2 


a  Vectors 

yS  vectors 

al 

Q'2 

q'3 

aA 

a5 

a6 

al 

>01 

0 

0 

0 

0 

0 

0 

0 

>02 

0 

0 

0 

0 

0 

0 

0 

>03 

0 

0 

0 

0 

0 

0 

0 

>04 

0 

0 

0 

0 

0 

0 

0 

>05 

0 

0 

0 

0 

0 

0 

0 

>06 

0 

0 

0 

0 

0 

0 

0 

>07 

0 

0 

0 

0 

0 

0 

0 
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Table  B.2:  Type  II  Error  For  P3 


Table  B.4:  Type  II  Error  For  P5 

a  Vectors 

P  vectors 

a\  a2  a3  Q'4  aS  a6 

OrV 

ySl 

0  0  0  0  0  0 

0 

yS2 

0  0  0  0  0  0 

0 

yS3 

0  0  0  0  0  0 

0 

yS4 

0  0  0  0  0  0 

0 

yS5 

0  0  0  0  0  0 

0 

>66 

0  0  0  0  0  0 

0 

yS7 

0  0  0  0  0  0 

0.070 

Table  B.5:  Type  II  Error  For  Pg 

a  Vectors 

P  vectors 

ckI  q;2  0-3  Q'4  aS  ab 

Qr7 

>01 

0  0  0  0  0  0 

0 

yS2 

0  0  0  0  0  0 

0 

yS3 

0  0  0  0  0  0 

0 

yS4 

0  0  0  0  0  0 

0 

>05 

0  0  0  0  0  0 

0 

>06 

0  0  0  0  0  0 

0 

0  0  0  0  0  0 

0.063 

Table  B.6:  Type  II  Error  For  Py 


a  Vectors 

/3  vectors 

a\ 

Q'2 

a3  aA 

OrS 

a6 

al 

0.006 

0 

0  0 

0 

0 

0 

pi 

0 

0 

0.005  0 

0 

0 

0 

/33 

0 

0.007 

0  0 

0 

0 

0 

0 

0 

0  0 

0 

0 

0 

/35 

0 

0 

0  0 

0 

0 

0 

/36 

0 

0 

0.007  0 

0 

0 

0 

pn 

0 

0 

0.003  0 

0 

0 

0.190 

Table  B.7:  Type  II  Error  For  Pg 
a  Vectors 


P  vectors 

a\ 

al 

Qr3 

aA 

a5 

a6 

al 

ySl 

0.532 

0 

0 

0 

0 

0 

0 

Pl 

0 

0 

0.453 

0 

0 

0 

0 

P3 

0 

0.589 

0 

0.004 

0 

0 

0.346 

PA 

0 

0 

0.042 

0 

0 

0 

0 

P5 

0 

0 

0 

0 

0 

0 

0 

P6 

0 

0 

0.020 

0 

0 

0 

0 

B7 

0.0004 

0 

0.100 

0 

0 

0 

0.070 
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